{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 频率分析 (2)：从 Gaussian 的 fchk 输出给出模拟的分子红外光谱"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这篇文档将会接着 [上一篇频率分析文档](freq_1.ipynb)，从 Gaussian 输出的 fchk 文件，绘制红外光谱。\n",
    "\n",
    "但为了文档的叙述便利，这次换使用具有高对称性的甲烷分子，而非先前的假想分子。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning\">\n",
    "\n",
    "**注意**\n",
    "\n",
    "不处于能量最低结构的分子一般来说不适合用作频率分析。此时绘制的分子光谱图从理论上是与不可能与实验相符的。\n",
    "\n",
    "这份文档尽管使用的是没有虚频的分子；但对于有虚频的分子，必须先进行结构进行优化，再作红外光谱分析。\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们先了解一下我们所讨论的分子与其 Gaussian 的输出："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAR0AAAEOCAYAAAC0KMNOAAAAAXNSR0IArs4c6QAAAARnQU1BAACxjwv8YQUAAAAJcEhZcwAAEnQAABJ0Ad5mH3gAAHqrSURBVHhe7d2Fu3PtVa79/C2f7U+2IntjG9gbK1CoUNritoEXd3eH4u7ubsXdKRSpG3Xa0hbq9tbyjTMr18rIyDVm7shc8rxPjuM82vdZWUnWSuYv933PObMWi//0IcvpPnh//1F90GAfONZ/oA84oIfv79/Tww7oofv7/9T7D/aQ/f2/6v0O6MHT/T/qQQf0wOn+79wDBnvf6f5d7n0O6P7T/V/qvQ/ovab7P9V7HtD9pvs/cu8x2LtP97/n3u2A3rXvf8u9y3iL//ShAYfDpithM9W1I1Sw6ZoNIYdNV8Km62CEEjZdsyLksOlaYzPVwQglbLruIrTOYdO1Bsc1itDiPwc6udkQcth0FWy67iIUOXBcCZuugxFK2HTNipDDpith03UwQgmbrlkRcth0rbHpOgqhNTZdHUI76LgsNi6DTZfFpqtg02Wx6SrYdB2EUMGmazaEHDZdCZuuORCiW4WQw6YrYdN1MEIJm65ZEXLYdK2x6RJAi//8YQELGWy6LDgug02XxcZlsOmy2HQVbLosNl0Fm65zrwfRbAg5bLoKNl0HIVSw6bqLUOSw6UrguM6J0AadmsHG5rDpMth0WXBcBpsui02XAcdlsekq2HRdO0IFm6ksNl0Fmy6LTVfBpuvc60E0G0IOm66ETdfBCCVsuk5ByINTc9h0OWy6DDY2h02XwabLYuMy2HRZbLoKNl0Wm66CzVQWHJfBpsti02XAcVlsugo2XRabroJN17nXg2g2hBw2XQmbrkMQWvyXD19ucuC4HDZdDhuXw6bLYdNlsLE5bLoMNl0WG5fBpsti02WwsTlsugw2XRYbl8Gmy2LTZcBxWWy6CjZd92WEyIFD2+jkApe7CA1msOmy4LgMNm0OG5fDpsth02WwsTlsugw2XRYbl8Gmy2LTVbDpsth0FWy6zr0epCw2XQYc1350aoHL2REKXK4VIYdNl8Omy2Bjc9h0OWy6HDguh02Xw6bLgeNy2HQZbGwOmy6DTZfFxmWw6TLYPOc5z1k+//nPX774xS9evuxlL1u+6lWvWt57773L17/+9Vvp397whjes/j9dIOSw6TLYdFlsXAab2uK/fERgohw2XYHLrUAocLlWhBw2XYHLMEQOG1fgcq0IOWy6HDZdDhyXw6bLYNNlwXEZbFJ///d/v3z605++guZFL3rR8iUvecnyla985fK1r33t8o1vfOPyTW96U9ub3/zmZXfhaxejIYdNl8HG5rDp2otOzWHTFbgMIeSw6Qpc7iIUaOzLYdMVuJwdocDlWhFy2HQ5bLoMNjaHTdf9l4961KOWj33sY5f/9E//tPznf/7nS2xe8YpXLF/zmtesRjD1AiK50csGH+Ww6XLguBw2XaDzFoGLsvDkHDauwOVaEXLYdAUuZ0cocFEWnpzDpiuAuYvQQA6bLoeNy2HT5bC56DGPeczyyU9+8vKZz3zm8gUveMHyX//1Xy+nUa973etWo5u5LoyOtgE6BCGHTZfDJrV4i48McHKBy61AKHCZE6G9EDlsugKXsyMUuFwrQg6brsDl7AgFLrcIIUY3j3vc45ZPfepTl89+9rMtOKzPXNVlF581QFeB0C46tcBlCCGHTVfgchehiyw8ymHTFbicHaHA5VoRcth0BS5nRyhwORGhP/iDP1j+7d/+7fLxj3/8jQFHFw9PLnCZA6HFW3xUAmZfgcu1IuSw6QpchiFy2LgCl7sIXWThyTlsugKYa0PIYdMVuAxC9Gu/9mvLP/7jP17+zd/8zSU4z3rWs5bPe97zVnumXvrSl14uGAPOIes057x4bLoCl3MgdIFOzYHjClzOjlDgcq0IOWy6ApezIxS4XCtCDpuuwOXsCAUutxyhRz7ykcs/+qM/ulwwZg3nGc94RrtozFrLdV82u+0dNl2ByxBCFZ23DGTI4kMOm67AZW6EJiFy2HQFMNeGkMOmK3A5O0KBy7Ui5LDpClzOjlDgkrP4KIdN1wOWv/qrv7r8wz/8w+Vf/uVfLv/hH/5h+cQnPnG1l6qbVs25cHzoJR8vNCtCi7f8Xxfo1CxA5LDpClyGIHLYdAUuZ0cocLmL0CaLj3LYdAUw14aQw6YrcDkDQr/4i7+4/L3f+73ln//5ny8f/ehHrxaOn/KUp6z2VGmUw7Qqj3Kua1rVXTw8OYeNK3DpELpApxa4nB2hwMV1ZQg5bLoCl7MjFLhcK0IOm67A5ewIBS5zIzQJkcOmK3A5EKGf/dmfXf7Wb/3W8k/+5E+Wf/3Xf70a5TzpSU/aO8o59Libq7gs/l3goiw8ymHTFdgoj04tcHHdRajksOkKXM6OUOByrQg5bLoCl2GIHDauwOUaEPrJn/zJ1cIxe6r+4i/+YmuUw1rOc5/73OW//Mu/LP/t3/5tZ5RzE9Hhsn3kdABzToQWb/nRAQg5bLoCF9fRAFHg4joJocDFdRchk8OmK4C5NoQcNl2By9kRWsOz7kd/9EdX06rf/u3fvhzlcGrDE57whOXTnva0yz1W9RQH9ljpVIbbAY8B6FiEFm8V4NAlPsph4wpcuq4EIYdNV+DiOgmhwMV1ZQg5bLoCl7MjFLi47gMI/eAP/uDyZ37mZ1ajnN///d9f/tmf/dlqF/k//uM/bk2tdAKndpPruJybjg4XD08ucDkUoUt0ancRSjlsugIX10kIBS6uK0PIYdMVuLhOQihwcV0ZQh6dH//xH1+Ncn7zN3/zco8VUyvtJudETs4a7w4G1AmdNxUcXTw2XYHLPoQWb/UxgYwKXLqOAogCl66jEQpcuu4iVHLYdAUwritByGHTFbi4TkIocHE1CP3wD//w8qd/+qeXv/Irv7KaWnEg4F/91V8t/+7v/s6u52R07r333luGTj2HzGHTFcjUttHJrbHpmh0hh01X4NJ1NEKBS9eVIOSw6QpcXCchFLh03ccR+oEf+IHVKOfnfu7nlr/+67++/N3f/d3ln/7pn16u53AEcofOq1/96p3p1W24TJ9J77Bx7UWntsam6y5CKYeNK3DpmgMhOgogCly6jkYocOm6EoQcNl0beH7oh35o+RM/8ROrqdVv/MZvXK7ncAQyu8p12gPH52R0Xv7yl++gc5suK3hqxyC0eOuPDTSUw6ZrjU3XUQitsemaHSGHTVfg0jUHQnQrEXLYdAUuXUcjFLh0HYHQ937v967QYWr1S7/0S6v1HHaVc0Agi8igw54rh04e6dyU0x8OuWyfzBq4uEYAWqFTmxUhh03XGpuuoxBaY9M1O0IOm641Nl1HIbTGpusuQimPzo/8yI+s0PnlX/7l1QGBGR32XE2hkxeSbxs6XLbh2QMQOYQsOrW7CK1z2LjW2Ew1K0IOm641Nl1HIbTGpmt2hBw2XYFLV0Hoe77ne1bosIjMrnIWkUGHPVdupNMtJNfjdG7TxaNTC1y6gGfx1vcELLmETdfBCCVsumZFyGHTtcam6y5C6xw2XWtsuo5CaI1N1wwIgc73fd/3rUY6nPbAyZ3suaro1IVkdpnXj7NgeqXTIG7TZXM2vcOmq8Czi04tYdN17lEQHYxQwqZrDoToIIAoYdN1MEIJm667CK1z2HRtwGGU8/3f//2ro5BH0NFxOvlD1/NpENplftsu7uM8PDZNi7f+uIDDYdNVwHHdFoTonFMxNRtCDpuuhE3XwQglbLpmRchh07XGputAhIQORyH/2I/92CU6eXrF3it2mXOcjg4O5IjkfBqE9mDduz5W51ZOsYRObRShxX8NdHIHIVSw6ZoNIYdNV8Gm6z6PkMOmaw3NVAcjlLDpmhWhHhymVuy50jE6WtNhIZld5hyno4MDdRoE517xsRYvfOELVyd85sVk7cECHQfPYrHYyV1Gr3fOy8XZ9WtsujqEdtCpWWy6CjZd514PUhabroJNl8Wmq2DTde71ILqLUOSw6UrgdK0RAh1iaqVjdBjpgA67zDlOh4MDdUQyp0Hw4V2c8Kk9WHVdJ0+xbi86tTU2XQJo8V8/PnAhA47LYtNVsOkyCNVzVhiS8iQxJOV/p/7/4QgZbLosNl0Fmy6LTVfBpuvc60E0G0IOm66ETdfBCBVsSkKH0Q7osKbD3it2mYMOH9zFGeagw7lXebc5i8luXUdTrDsLnVzCprZBJ2ew6bLYdBVsUjxR+S8d5jNy9cQwFFX5v90Tpgtfu4tQZLHpKth0nXs9iOZAiCw2XT06jHbYeyV0OCKZ0yA494rP0nGLyVrXqVMs7cXqXsOjmIxe75yXzecJOXBce9GpGWy6LDauj738W0BTf+nQPRn8mxq9cF3wugTorRw2XQabLouNy2DTZbHpMuC4LDZdBZuuOxAhgSN0SAcHchoE5179zu/8zuqD2OtiMus63RRLr/F71wvK7rU8isno9c55WW1L+bOFDgFo8d8+YXmRw6bLYNPmwLlnNbLp/tKhDp6a87LBZ06EHDZdBpsui43LYNNlsekq2HRZbLoSNPuy2LgMNl0GHKrosJjMwYE/9VM/tfyFX/iF1WfpsNs8LyYzxeIvebKu46ZYjHbYi5VHOxrN58soJqPXO/dlC53aFEIbdGqByzBEDhvXPZejG54MVvavA5x82cVHOWy6HDguh02XwcbmsOky2HRZbFwGmzaHTZfBxuaw6TLYdBV0NMrR9IrF5J//+Z9f7TbXYjLrOnymDn9cT1Ms3ljzXiwdnZxHO925WA6T0a7icvG5QgWbroyQB8cVuAwh5LD5uNU7AU/GTQJHF49OzWHjcth0OWy6HDguh02XwabLguNy2HQ5bLocOC6HTZfBplRHOnysBcfqsNucPVh5XadOsdiLVadY3WjnXjPNcpiMNveFx5k/0GxTwca1+G+fGJgoh01XADOAEO8EPBkssAEOi2s8Ae4vHXZrOFd18di4ApdrRchh0+Ww6TLY2AIXsujUHDZdDhtXhWYqh03XNDoZHg4QrOs6HCSYp1jai5UXlDXayWs77MnS6z+Pdhwmo8196dGp7UWn5rBxBTAGIea7PBHIzxw3/6VDVvLzovG9IX0dXl7HZXvXvQPHFbjccQg5bLoCl7MjFLhcK0IeHdIUi3UdTbE4SFBTLO3FcgvKWtvJe7J03M69abSz2qgNJqPNfVk9vssPNnPYdIHO2wQuysKjHDZdH786cIopFXNcfvk6WEofVJ2nVdotflMu2/AcChAFLmdHKHCZC6EhiBw2XYHLrUAocNmDUIcOi8ms63Bkcp1i5aOTebPNox2t7WhPlk6NyNMstgXegEcxGb3eOS+rx5c+2Gw7h01q8TafFODkApgTEWJVnwOm+OVrHYdfejfKue5plbt4eHIOm67AZQghh01X4HJ2hAKXa0XIYdMVuFwBQg4dretwkKCmWHkvlhaUdcxOHu1oT5aO28mLyvqAL0Y7q13SBhO3nbjrzX2ZRicXyFSEdtHJBS4HIsRBU6zmc7AU2vNLz3/PmVEOwtdRjoaUN+myeOvARVl4cg6brgBmLoT2QuSw6Qpczo5Q4HKtCDlsujw6edc5f2iPKRYj+7qgnNd28p4st6icp1nam+UwuQno8BhWKG59wJkDx7UXnVrgMoEQw03mt8xt+aUzysl/joNfNItoWrnXoeACR92ky/aR04HL2REKXHIWH+Ww6Qpg7rMIBS5nQqhDh3QeFm+0fHQpC8putKM9WTpuJy8qM83SQbEVHoeJ2z7c9ea88BhWKG6hk3PYpBZv88kBSM5h07XGJ2JhjWkVC2r6cxz8shGeXzS6d1Mrhmo3FR0u2/DkApchhBw2XYHL2REKXK4VIYdNV+BygxCaQkdTrLygXNd2dD7W1DRLe7Py+g7bhsPkutDJ2yejHAYMi38fwJCFJ1fReduAhnbwUQ6b7ZjjMrfll661HPZYMazUAjILaPpckbybnPnrTUeHi0enFsAMQeSw6QpchhBy2HQFLmdHKHC5VoQcNl2ByyBCDh2lBWUds5NHO3lPFttCXlTO06y6vpPhcZjUbYXc9c590X1x/wwUeIz5UxUvARpBaPG2n3KBTs0CRLvo8EtnlIPyDC35RTO14peM7HU9B3S0nqNFM/0yb/Ll4hQOh40rcLmL0CaLj3LYuAKXK0bIYaP0+Tp5tJPXdjhchI+8qNMs7c3Skcpa38kLy6uN2mCibUVdBToCh/tme2WwwOPcRqcWuHQIXaBTC1yGEPqky/NQ9Mn4TK34JTOkzOs5iK7dhHkR+fahk3PYdAUuZ0cocLlWhBw2XYHLLUTIYaOYYjHK1wd7sbajPVlsC4z686KyTo9gmpXXd6bg0VRLu9O1vWibuYoL2yX3xX3zGHhMTAk9Nl2BjfLo1AIYgxDDS310I79odGdIyS9YB0ZpV7nQ4Rda0bkN4OiyOoG1OYnVY9MVuAwh5LDpClzOjlDgcq0IOWy6ApczI+SwybENsLbDmiZ7snTcziMf+Ug7zcp7s5gFTMGjqVbdqyV4rmqbYfskwOFxsP2yPddPWPTYmBZv+6mBCDlsXBt4dA4KUyt+ycg++udVQScf+n1b0OFyCU/uLkImh40rcLlWhDw4ymGjmGLpuB2dGsH2kBeV8zQrr+8IHtY8HTxaXNZeLZ0uoRHPVWw3eZQDejwOtmEe9+I/VHRqBhxavJ3QyTlstuMXzi+aqdXUn1fN6LiRjtC5TZftk1kDF9dJCAUuritDyGHTFbjcVoRaiMbRYYoleHSUshaV6zSLbUPrOxWebsQDPDqOhz2+mm7lUY9GPmxHgkidehE63JdGOTy2zYecBS45iw9VdGoDCPEL13pO/vOq+9DJazq3daTDZRueXODiOhogClxcJyEUuLiuDCGHTVfgcnaEAheXQahCUwMeTbNYVM7TLPZmVXjYja6F5Sl4tDtdBxBqnSePegiAGP1UhPLaTy6jtC++X+Bw3zwWtuf8yYrbBTD7EFq83acFNBS4dBWE9ItmAY3hpBaRp9DJC8l5lzk/ED8YP+BtunhwXAGM6z6BUODimgMhOhogCmBcA+gQox3gYZvQNEt7s/L6DjtaWFgGHhaWp+DR7nS2GzZ2rfPkUY9GPsInA6QEUcWoS9fle7k9gAMc7h8M2UHkwakFMA6hDTq1NTgm/ZIZ6bBizy+2opPXdPgl8gtkWKbjdPLBgfygtw0dLpsz6h02rsCl60oQcth0BS6u+yhCFRlXnmZpb1Ze3+HjL9hGMjzdiEfH8ejI5Trd0qiH7UgjH3IA5TJGuXwdfS+3BWogBzg8Brbl/MmKHpuuNTweHNcGHM1jWbGfQoe9V9plruN0QAet+WXxS+KHRFWGfbft4j7Ow2PTFbh0HY1Q4NJ1JQg5bLoCl66jEQpcuk5AqALTJXj06YKs77DmOQqPdqfzZg08zBLYdjTdYtQzhQ9ICB8BpIRJLV9H38ftcJsCh/tl9AWMi/8Y2NAWPsphU1q8/acHJsphsx2/WNDR/NWhw94r/f0f5NbBgTryUude8YPxQ2vYd+umWEKndhehlMOmK3DpOgogCly6DkSoAtOl7UO70UfhYZthdzrw6MjlvGeL2YJGPVrryfho2iV8BFBGqKavE9fne4UNgwNuH3DYfoGQKeElOrVRhFbo1CYQ0ihHuwj5hbKm0/3RMdDRaRD8wrrd5kLHwTN6xOXo9c55uTi7fo1N1+wIOWy6ApeuoxEKXKa6QxBywHRpG8kLy9qjJXi0xsMbtXan6zwttp063WLU4/DRyEdrPgKojoBq/LuQIa4vbLgtbpf74P64fzDMn6q4KXBxdQhZdGprhOovVPNW0NEucw6EAh1+eQwXdeSlTufPi8n8oGjLFOt2o1NbY9N1FEJrbLrmQIgOBojW2HQdhdAam67ZEToMHcp7tNjLK3jyiIdtJu9OB568wKzpFqMeLTJrrSfjo2mXRj+gwQgoI6TARenfuB7X5/u4DTBjZsL9AA4AMqDw6NQCly7gWbz9ZwQsuQJOSr/MPL1ilZ5fJL9EdgsiN788fnE68tItJmuKxbBOU6w7B53aGpuuuwitc9h0rbHpOgqhNTZdgU9GZV9sJ25Xep5qsVdLu9PZqDlyuU63NOrRWo+mXOCjkY/WfASQRkBCSBC5+BrX4/p8L7fDbXL7QMd9g2H+cDOPTVeBZxedXI8Oox1W6YUOuwXzBxjxS9N5Jjqr1k2xGN5piqVjCOplFJPR653zwggtf4jZWGtsug4GiNbYTHUwQgmbrvsgQhmWfeXlCMGT13jYblgP1ZHLzBTqdItRT17rYcolfPK0S6MfAGI7ywgptj2lf+M6XJfv4zbAhttluwU8BhA81stPVjSfsOixca3Q+cw1MNPVXyRxTALDRnYLuk/E17oOQvMLqlMshnRMsfJox+3FGsVk9HrnvuTPFdp8qJnDxpWw6ZoVIYdNV8Km62CEEjZdsyLksOk6HB3q4GG7YScM2w5ronmBuY56hI+mXOCTp10a/QCQRkBCSBC5+Jqg4XvBBtC4fe6L++bxXHy+0BqbrlGEFu8Q6OQMOFR/iUKHX2D+RHx+aVpM1sltOpWfXwo/IKoylGOKxWiHRS1GO9p9Xkc7o5iMXu/cly10anMgRLcKIYdNV8Km62CEEjZdgwjpw+ry9jCS4Mm701kPZdthTVQLzJpu8eadF5l5E8/45GmXRj8ZIPAQQgpYhIviOlyX7+P7gQzUmKUw0mLkVT/YbNMam64OoB10ag06zFV1EBTzVMRGa63roPTUFKsb7WhPVr44TEa7igsf8XFRwmYqi01XwabrIIQKNl13EdoChzdZjerz9jCathtGPMDDtsPyRF5g1nQrLzLnKVfFh+2LEUkGCDiEEG/2BCo1/l3Q8H1sp2DGbXM/3C+Py4NTS9h0CaHFO3yWx6bU/fLYbY7Ybl2nTrH4QdGVIR3zyDzaYW1He7LqNMthMtrcF0ZlG3RyBpsui01XwabLYtNVsOk693oQzYaQw6YrYdPEa51RCeBo/SVvE6Np29GRy2w/ebpVRz3cF7OHjI+mXdrTlQECDSEEIkAkjHL6d67D9fletlVA4/a5Px7H7oebVWy6Eja1C3Rq0+jkX572YPEL40Hyi8pTrLwXa2q0oz1Zmmbl0Y7DZLS5Lz06NYNNl8Wmy4Djsth0FWy6LDZdBZuuc68H0YkI8TpnNMLrWwgAABtn3S5GY6rF7TJ6YhtiuqVRD4vMWuvJ+DDy0bSL7Ys3do1+AAgwhBCAAJEwqvHvXIfr8n3cBj8PPxezFe4/f5jZTgcBRHvRqe2io18cWmtXoJti8cvgB3SjHa3tsCeLYwkY7bhd6A6T0ea+rB7f5WcMOWy6DDZdFhuXwabLYtNlwHFZbLoKNl3XjBCvcV7bTH200CtwNNVx28ZoWmDWdIttSWs9zB4yPtrLBT519ANAPCYQ4s1eEAkj0n8TX+d6fA/fz/bK7XIfqynV5Ud7FGy6DkFo8d8/e3mRw2aT+4UhtdZ1GB5qipX3YuUF5TraYfWc0U4+bodF5XzAINOsUUxGr3fOy+rxpQ82285h02WwsTlsugw2XRYbl8Gmy2LTVbDpsth0FWy6JhDidT0FDiMLRgxu+xhNC8yMeoQPe4W5b4cPb+x59COAgIPHxyiIbU8Y5fg3vs71uD7bKYMEsOH2ub/8WUK7FWy6phDaoFMbQ0er8dp1zhCUH4BfhkY7PDl1tJP3ZHHMQF5U1t4sLSo7TBhh1Iu73twX0GGq6NGpOWy6HDguh02XwcbmsOky2HRZbFwGmy6LTVfBpiuBw8auJQNez2y4GkUIHE1Z3DYyEksVGR7eyJlyCR+NfJh2seYDghUgUAQhABFEgFLj3/k61+V7+H5uC9iY3l18nIfDpqtg05UR8uDUPDr6ZfGLYk7KMFSnRNTRTl7b0flY7L7Li8qaZtW9WQ6Tm4QOj3vzAWcOmy6Hjcth0+Ww6XLguBw2XQabLguOy2DTZbFxGWxSbOiAk99AAYfXsgOHeH277WQktiVtT3ozzyMfpl2gkEc/bGs8Rt7ogYNZhiDicRMgKf0b1+H64MXtcJssaNeP87jMYtNVsHEt/vvnBCrKgXPR1C+KOSm/GH4paJxHOxqO8kTxJLFSzuo5u+u6aVbemwU8DpObgs5qJJY+UXE3h40rcLlWhBw2XQ6bLoONzWHTZbBpc+C4NuDwOmaD1mt4FBzltpWRMjwKgLTmAwxacAZFAcQMg+0uQ6SARf+fr3Edrsv3cRtst+Bmsemy2LgMOLSNTm4cHU2x8oJyHu3kPVk8Qeymq9OsvDerwuMwuS50uF8FOEwDF28XuJBFp+bAcQUuZ0cocLlWhBw2XQ6bLoeNK3DZgxDTGTZORgsCR6P1EXDyXiG3zexL21StTr0EEI+XBJEwyunfdT2+D2y4LW67nk1/kcGmy4LjEjrvGMCQhWeT+wUpLSjrmJ082tGeLM2F86Kypll5b1Ze32FhGXgcJmzwGQBy1zv3RffFtIqpH1PB/KmKFwUuQwg5bLoCl/s6QsMQOXBcAU1CiI2RkQDg8LplcdaBk3dDd+Cwjsn3ue3l2Bj5ZHx4oychRPn/6791PeL7GDmxzS7ewmHTZbCxOWxKi3f83EBHBTANQu6XoPhlaAio0Q7vFtqTpVV/dudpUVnTLO3N0vqOg8dhogMI1VWiw/1x/4xyeIy76OQCl2tFyGHTFbicHaHA5VoRctjsxjs/b5aM0DtwwKQDh/8WOFyX7+ONlte922aOjZGJpl0AxHYHQvyv4t/5uqZoudXIJp28up3DxuWw6dqLTm0NUOR+AYofhB+SHx5dGcIxz2Q+ybsGK+X5SeRJ4cnT3iyt7wAP6zssLLsRD1MtHTyoPVskdOa+ZHB4DCx2swB+8QFnDhxX4HJ2hAKXa0XIYdMVuAxD5LBxBS4nIMRoAHB4k8zg8CYJOLxRHgIOI3qBw+1xu2wTbts5NuGjhAz/nqvf57HpcuC4HDZde9HZrv4AOX5glGX4xjCVuSRDVVbJ64Kcplk8WVrfAR7Wd7SwrD1agkdrPDpqGXTY+IXOVVzytAoAwZCRWf10RY9NV+ByF6GLLDw5B44rcBlEiNE5b5DAwJuj1iC785wcOHwtg8OOE8DhNc/IiW0hL1ST24aOqSLjoKF8Dtl2DhtX4HIuhBbv+HkBSs6DQ+6HUfywaAs8vHOwaMUKOYvKdZqlvVk8UTxxwMP6jhaWMzwa8WhxOR/Hk0c8gDDnaEfgcF/sUeMx8JgYoS3evqJTc9h0BS5DCDlsugKXsyMUuOQsPsph0xW4DCHksOkKYAxCjMwBhzfGc4PDax1w2AYqODm3LZ2rzakcgQtZdHIOm67A5ViEFu8U0NAOPmoMHWK0o4Uu5shaVM7TrPzEAg/DUp5ALSxXePIaDxs5x/HoyOV8rpZGPnn0Q8Lo3CDxGHhcmw85C1xyFh/lsHEFLnMjNAmRw6YrcDk7QoHLTAjx5ggI5waH2wAcXu+89vOesU0eIHLb1Wj5fLHpApi5EBqBaPFOn3+BTs0A5H7QnEY7wKNFZU2z8t4s5sx6gnlyO3jy4rKO4+FAPK3zaNSjtR5GIAIoT79ywsih1DV1yZ+suF3gchehTRYf5cBxBS5nQIjXKDs7eE0CTh6BHwIOr1vA4c1T4IAX4PCaH9lF7+Dp+6hN6VQNn8PGFbhcNUIX6NQCmQYhh01O8DDa0TRLe7O0vlMX6zI8eaqV92rpOJ68zqNRD2s9eeTjAMojoTwimmr04tGpBS5DCDlsugKXYYgcNq7A5VoRcth0BS4HIsRoXODUKX8Gh9ejwMnoOHB4DQscXt8sKewuWAcw14aQw6YrcBmGyGHjCmhyHp3cNj4OmhrwaJqlvVl5fYdDsllkm4In79XiOB4OIMzrPHnUo7WejI9GPwJICKmMUU5fP+bisXEFLteKkMOmK3A5O0KByzUgxOuSkTfT/Q4cIDkGHF7LPTi1wOW+jNDinb5gjctYDhmXRjzszdL6Dot2h8LDcTw6gJB1njzd0qiHtZ4OnwpQhSiXv37KZfMxrw4cV+DiuouQyWHjWsOTAhxef7z2ACfvTQWQUXAYCQkcXru8hrlNMNssVDtsugKXsyMUuNQsQOSw6QpchhAy2KjFO4NOzYOjHDIujXi0G13w8OR08PBk8qTyJHMcDwcQ1nUejXqm8MnTrgyQEKrpa7o+o5196zlTlw08hwJEgYvrJIQCF9eVIeSw6QpczowQrz3AqYdv1KONR8DhzVHgcDuMmnhNszid945t57DpClyGEHLYdAUuZ0cocDkGoRU6tT0IOWC6BA8Lyzp+R3u02KWY4WGYylA3w8ORyywwa7pVRz2s9WR8NO3Sni4BpBFQhoj0b8R1mLrxvfx/MGLt55TL4h0CHLoShBw2XYGL6ySEAhfXlSHkwfmUT/mU1ZriMeBQBYfvy+CwPsQb69aueWUBIodNVwAzBJHDpitwuS6EFu/8hQENBS5dBiEHTJd2pec9WoJHI568V0tzbL0YeOKZbuVRj9Z6HD5a82H0I4AYAWnxGYhcfI1RFOtGfD//DUZMvdjTdezlEp7cXYRKDhtX4OJqEPrkT/7k5Wd+5meu9iidGxz2fAEOO0zy3rFNgYvrJIQCF9eVIeSw6QpgHEIbdGqBS1fA43Dp0voO8PAECZ68xqO9WuwBYHFPLwyeaDfqYa2n4qNpFwABhwDSCEgQCaMc/8b9cPvcJrfD9/M1RkpMvVhoPmnK5fBRRyMUuHRdCUIOm67AxTUTQp/0SZ+0AufzPu/ztvaYVnCAZRQcXpcCh72xrFd6cFyBi+togChwcZ2EUODiOhdCHhzXLj4OmC7gySOevMYDPNqdruMmeGK1zlNHPaz1OHzymg+jHwHECIgEkYuvcZ+8sHjxMaLitvg+0GKUxJSLReZDdqe7y+YD7wOXrqMAosCl62iEApeuK0HIYdN1gc4nfuInLj/jMz5jBc4XfdEXLb/kS75kBxww6cDhvwUO1xU4vCECGG+UvHnWPWQem64AxnWnI7R45y8KQJTDxnU4OtTBw+505twMVfU5PFrnyaMerfXwQqn4gITWfDT6AQ2NgASRMFL6N74OfLygeHHyzsbtAxpfByZGQ0zDtNB82pRr+69t3EUocOk6EKFP+IRPWIHzuZ/7ucsv/MIvXIHz5V/+5ZevIY2eR8HhTa+Cw2Egbg/Z8QgFLl1XgpDDpitwcY0itPgfgU3uQIAcLlNleLQ7nSeQA7Y4ilPrPHXUo7UeXjj5RcOLhGmX1nw0+gEgwBBCjIKAyMXX9JEcjLqY/zPNY5TFfXC7XIdRkxaamXKdZ6HZ4UNrbLpmR8hh0xW4dF0xQg6cL/uyL1t+5Vd+5fKrvuqrTgKHN0LeFHmj3F6sDly6bhRCgUvXlSDUoVMbQMjhMlVe4wEeDiDkyGUtMGu6xcavUQ8IaMqV8dG0Sy8eRicApBGQRkGCSBiR/puvab1JH8uhI1Z5sfHuCGxcj5GTW2g+Zcrl0amtsek6CqE1Nl17EPqu7/qu1fOpI9AZvfImAuD8HhUbqeK/ea4PA4gCl641Ph//8R+//PRP//Tl53zO51hwvuZrvmb5tV/7tZPg8LUMDq+1DA5vTNvguAKXrqMAosCl67YhtPgfXxy4qIRNl0GoojJSHvHolIk83aqjnil88pAZgPIISAgJIhdf0+PSyaoAqA8i4355t+P2uf58C80Om641Nl0zIPQt3/Itq98Rz5s+xkTPG5Dw3LFR8vwBN/EmwnPJ75I3E55TBez8L19bnatksenahufjPu7jtsD54i/+YgvO133d1y0f8YhHDIPDCJvRNm+A/Dx5z9hYa2y6ZkfIYdMVuHQdjVBCR22jU0vYTBUIaaM9NI16eBHrnVLna/Ek57UeTbkyPrwoAKECxAtHCAmi/A6nFxnx9fyY2Kh4LDwONiA2DO6bNSbug++dWmg+fcqlP/3jsOlaY9N1MECfvtpgv/Ebv3H57d/+7avfi54njU7BhpELQOv5ysgIFkat4M1zyMZb49+J63F9vveQkdA999xjwfmKr/iKFThf/dVfvQWOquDwOuB1w2tI4PA64zGuUMx7x+gSIHLguNbYdB2F0BqbrtkRcth07UWnZsBJ5Q330MAnT7c0ROfFrClXxoeRD4t6vAuBQQWI9RiA4B2LeCHV+Heuw3X5nvqY2MB4LLyD8xi00Mz9AFteaGbKddaFZv7umDoYoYRN1wRCbKBg823f9m2rKVQd2eiNQSMbYaPRTIWG54s4HovnLcebCOV/47ps6NwGtzu1HvSxH/uxy0/7tE9bfvZnf/byC77gC4bBoa//+q/fCw6PhZ/pYs1ojU3XwQDRGpuuOxGhxf/8kuVWFpuu88OTRz3CRy9uvbB5UfOC5oWpF7IA4oUCQszDeeEQUzEXX+N6XJ/v4zbcY2Jj412dx8ALkPvl+rw4ebGec6GZx6XHxM+4hc/MCDEF+YZv+IbVNOo7vuM7Vj9/Xq/JbwaMbhgF8pwwCtDUqWKj50e48DsmRqvEc+bia1yP7+E2Vr+LrcXpT15+zMd8zA44X/qlX7oFDj9TBw7/y9dASeDwBsTrgt//Nji1hE3XrAg5bLrW2HQdhdAam64phHbQqVlsXOeBJ6/zsLHXF7veWXmha/gugPQi14s7v7BBQulFTXpR652V2+V+6uNi45taaGaYroVmplzHLDQLSEZr/H9Gbzze1dDe4aNOnIqxobLxMbr51m/91uV3fud3rn5mPQ+a9mZwutEN2OQ3gwpNfi7YsImfk/hd6v/ra1yP7+P7Kz6f+qmfuvysz/qsk8Ehrs90kt8/98t98fry4LgSNl0HI5Sw6ZoVIYdN1xqbrozQ4n9+qcemy4KzW91oR2NkQeDDhs4Lv64jaGjv3mmFkN5thVGNf+frbCh8H7fBxsR9TH2MR15o5nu4HzYKLTSzKM1u+rzQzJRr30IzL3SNzgRPxoeNwILjGkSI3cpspGyE3/RN32TB0XSqol/B0e+b320FJ2MjZPhZ9fN26TpCiNsQPgLn8z//81cH/50DHI7l4XZAjNETa0WXi9Vl1/z+EjZd91WELtCpGWxsHhxVN9jRMjy5vK4gfPTOy8agvSaCiI1CG4bSv3Edrsv3ARgbFhvZxUmt058hxGPh3Z8NcWqhmWOAmHLtW2hmI2JjYqPUBslGJ3zylGv17nv5RxALNlMVcNhQ2UDZACs4I1OqEXAcNvXny8C6dB0hxPdzAmcFBzAyOGCSsaEMjtDpwGH3e947ttNsCFVopkrYdB2MUMKm61SEPDo5h03XPPjktOaTp14AxIYBQoIISBQbi/4/X9P1gIbv57aAJJ9Jv6kf9TAa4Pu5XTZERk1sbGwgeaGZo5210MyUKy80gwgbLRssGysbah4NcFvaAPn//DvXyX999RCE2MPD8StsaGyAAoc1nO/+7u++xD3/bvldZXA0ssxTqgzO1OjGYcPamGJNpaavcV2dwKnTG04Bh+tmcNjzxe+HAww3i9UFm66DEErQTDUbQg6broTNVIcgtHiXL1te5MBxOWy6TofHVfERQMQoyKWvK97J+X5uZ/HODhufeyzcvxaaGQEwdWODY0PJC81MufJCM1gx2tIGzMarUQ/fXzdUNjr+l39jQ+b6Hh+1DQ4bKmsfbGhseKzhsGjMLvHRvVR5Kju1hjMFTsUGWBgh5rTXMf9bBUenN+wDR1Vw+L4KDudsbcBxFWy6bgtCZLHpKth0OWzUBp1c4DITQnWjPTY2djaSvO6j9N/8L1/nehopKW5D55Dt5sHJ1cfD/bChspGycYIIGyEbHBuPW2jm+ozQGEUwgtCGzEasUQ8bbocP/87XgcCjs4kNCnBY92CDYy/VN3/zN692i58yrdKaWQZHj7sDx2EjZEiHNNQcOACawQEVBw7tA4ez0i/2jjlsuhI0U1lsugo2XedeD1IWmy4Djms/OrXA5cwI1Y322IQIsAihDpmcx6ZrDB4Cu26hmd2yWmhmKsNjZVTB9dm42bDZqLVBszGzITNq0AZcRwv8N//ORn/x11h3wWH9gw2V41fYQNkAu2mVRjl5wX5qlDM1rToEHMHCtLSrnsB5CjjgCzj8Xrhd1okud8nXLDYug02XxaarYNPmsOky2HRZbFwGG9fiXb48QCGHTVfgciaE3IZ7TBUZBw3l0ze2c9h0TePDBsxoQQvNbKhsnMDBRsbBaIwydCwMIzKuz4aep1wa9WjKNTXq4d+4fcDa/Dnoz16d/MjiKBsYGykbno7F0bSKx3DKKIfH141yeJwVHQeOYAHmLoEzdXpDru6pyuBwGwKHPWEWmy4Ljstg02Wx6XLguBw2XQYbm8Omy4BDG3RqDpuuwGUIIQ8P5Y323G2OIQpcyKJTc9h0+VGPFprZePNCs0YcbPyMNrSewkbP9RlhMLpgZHHsQrMOnGM6wobKhjY1yuGx5rUc0MmjHKEztZYzMsoZBYfDD1znBIffDb8jdzKrxcbmsOky2LQ5bFwOmy6HTZcDx+Ww6RI67xrAkIUn57BxBS5nQIjchjxSPlBxusBlGCKHzW71sTDiApO80KwFXTYANggQ0PlNbtSjKZcWmjXlmhr1MF1g1y+jA22o3Jdby8lTq7zHyk2teBwOndG1HKHTgVORYSFe8d/nAIffPc8Bo0APTs1h0xW4DEPksHEFLteKkMOmy2FTWrzrVwQ6KnC5wQht+uJNade8z2HTFbgMIeTByWV4KC80c3KiDt1n2sMGwUbClAcIgAoEjl1oZkFUn5pXRzl1jxX3NTK1yugw8srrOZpadXusRkc5HTY1fpZDwOH6FRwWpXfPpHfguBw2XYHLXYSm0KkFLkMIOWy6Apc5EdoLkcOmK4A5AaEKD/ERDOyWZR2BFz5rCgz1WRhlA2Ekwijk2IVmDmrT+UiMcrSWw22zQdapFegwugKdPLVy6zmgN7Wek9EZHeVohJPRycDokwBqx4ADwBxQyML6Lji1wEVZeHIOm67AZQghh01X4HJ2hAKXWRCKPDauwOVaEXLYdAUuZ0cocDkSoYzOR3/0Ry9ZbwEHjjth1KOja9koAIKNZnShWVMuNnpAy8ey8K4OZoykuM06tdq3nuPQ4T4zOm495xh09oHDUd5K/zYCDr/PXXB0Rr3DpitwOTtCgcu1IuSw6QpchhEig41avOtXBiI5B44rcDk7QoFLzQJEDpuuwGUIIYdNV+ByAEIZHcW5PYx6WH9hOsSxIlqDYYNhQxpdaAYybkujHCGmqRUbo6ZW3XoO6NT1nHOi48ARNh04GZtchqeCQwIHdAUOv5f8ER67OWy6Apf7NELksOnK6LxbRSfnsOkKXIYhcuC4Apc7DKEKD1gwQmEdBjB4J64LzWCRF5rdlEsjp4xXnlqxUdb1HG6PKdzIIvIh6LhFZI1y3AgnYzMCjsrwdODwO9D5VDzOxTsELjmLj3LYdAUuQwg5bLoCl7MjFLgMI0QOmy6HjWmFTs0CRA6brsBlCCGHTVfgcnaEApdrQCjDQ3wYFef8MOrpFpqZGrmFZkZMfK/WiUBLBwMKnbyec9Xo5FFOHeFUaBSocBQ3Hzlbscnp+vvAAXYeH495+6M9ApezIxS4XCtCDpuugGU2hCKPzlcFNLnAxXUSQoGL68oQcth0BS5nR2gNT0GowkN1oZlRy9RCcx0paWrFSAl0AIuN0aGjPVegoz1Xx6LTLSTnqZVGOd1USn9MkaO2BQ4nzeojRbv0/dxWBofRHr8/QAY/MARH1sO24ckFLvdphChwmROhXXRqgYvrLkImB45rG6AKT15ozlOuvNAMICw0M0KqUyve3fMickWH7zsEnam9V1PouKkVo5wOGz4ORH+pI3+o/qHwgBr3A3DcJ/cNfDweHhuP2YPjClyGEHLYdAUuZ0cocBlGiBw2XYHLMEJksMkt3u2rAxGHTVfg4joaIApcXCchFLi4rgwhh02XH/Vo2tQtNHOdPC3L6AAU6GgRWXuuOnSmplcOnX3H6Qgd7bXKUyuBAyQZG87E54RY4uRYwNHIB4S4vhqFB3S4bx5DHu2A5ubzhRw2rsDlWhFy2EwVuMyCEDlsuio67w46NYeNK3DpuhKEHDZdgYvrJIQCF9cJCFV4KE+fdGwPIxm+BjraA6avAVPdczWKztRxOt3BgYwcQIeNWeiwcTO9qus5YAAMQCJs+LAzoOGjP/Snn/lLG3ydjwUBHkY7FZ6Kj8Cpox2hw6gLDHmMPO78oWbbOXBcgcu1IkQOm67AZRghcth0OWyaVujU7iJUcuC4AhfXEQhVeDSi0bE9+nftds+LyB06U9Mr0Jk6ODAfkcy6zshpEFrTyes54AAkeWTDpysCDZ+wyMe78imL4MPX9QcTBY/DRuBMoaMpFo+Jx8fj5WdZ/PdARp0EEAUuZ0coYLmVCEUOHFq8+9cENBS4dB0FEAUuXUcjFLh0XQlCDpuuwMU1iFBGpyLD/2ftJ6OjtR/tLu/Q6fZedehoMRl0RvZgVXQYZYAOWAAIH+2hkQ3YvOhFL1p9sBmfJc1HuvLJiuADRvxhQ414NM3qwBE6gJPRqes6W1Os1QecBS6uK0PIYTNV4DILQhS4zI3QBp3cGpuu2RFy2HQFLl1HIxS4dF0JQvvxoYxOXnB26HS7zLsjkkGn24NVF5Pduk4d6bDxgweI6O+ECRs+RZG/msGf7OEvZ/AZ0vw3X+ejXlnf6dCp4IygwxSLxwmW/Fz1ExYtQHRlCJHDpitwGUaIHDZdgcswQmSwqXl0amtsuu4ilHLYuAKXqQxCFZy6uzyjU6dXQkcHB+YjknXu1cgerNF1nbqmAwDgkT+yVX+YkA+r5y9l8Gd6+GsZxL8xzeJTFhkZ1SlWBWcfOnVdR1Msfo5ddGqBi+vGIkSByywIkcOmy6LztQEHOWy61th0HYXQGpuu2RFy2HQFLl0zI5TR0e7yDp26y7yeBlFP+Kzo7FvX2TfFEjps8GDAOg7TKf09MEY2YKO/hCpwiNEOUy3WefKistDpwBE6gNOhw+PSFIvHfvEJiw6brsDFdRJAFLgMI0QOm67AZRghcth0BS6HQLR4jwAnNytCDpuuNTZdRyG0xqZrDoToYICooLMuj3KYWlV08t6rjE4+DaKis29dZ2qKNbXrnI1b6znAwIiFUQ5/A0x/DUPIVHQIlBgRMTp6ylOesjXFmgJnCp28rqMpFj9P/YhXj01XAOPaQYgcNl0By2wIUeAyC0JksFE76NTuIrTOYeNK2HSdgFAe5YCO9mhNoZPPvdJZ5tqDNbWYnKdYbtf5vr1YGR1GJ+yNYp2GqRNTKdARPLpkdBgFMQVjdKQ9WUyxtK7TgePQYYqnkU5Fh59j9dnSuRuJEAUucyM0BJGDZqotdL4ucFEJm66DEUrYdM2KkMOma43NVNeM0BQ6Ok5HRyTnc69AJx+r4xaTu3UdN8Vye7HqaIe1E0YWbPhMizgeh5EO06aKjuDJ6PA11ny0oMwUK6/rjIKjXeZ1eiV0Vus6+lD7ritDiBw2XYHLbBAFLsMIkcPGtI1OLWHTde5REB2MUMJmqtkQcth0JWy6GoS0niN46vSqQ6fuNs+LyYdOsTTa0RSrG+1kdNjwGZkwRWKqxAKxplcVnXphbxbXZ/c507O8rjOFTQYnT60qOiDJ4wfVzQfabz7Yvu0khChwmQUhClxmQYgCl1MR8th0FXBcdyxCBZuu2RDaRocY7bgTRDnhk5NE825zt5jMFKtDx+3Fmhrt1LUdNujFYtHGCEajnTy6qReuw8iIXewsRGt6xdSqYlPBqaMcreeAYUYHMFe7znfQqQUuXTsIkcOmK3BxWYTIYdMVuAwjRA6brsBlGCECnfs9YnnZe5DDxlWw6ZoDIbLYdBVsum44QhmcvJgsdHTuVUbHLSbX0yGmplhTC8rd2o5Dpqsb7fBv7MFiRMRiMtMrjXIYNTHKqdDsA0frOdplntHh51m84+dGDpuugKXrRiNEgcssCJHDJrWFTs1i01Ww6Tr3epCy2HQVbLosNl0Fm64T1oM6dPIJn3y0Rd1tXtd18kGCI1OsQ0Y7DpZ9AQwjGhaOQYaFZvZwsVudAwdZB+KgQkY42l1esRE0wkbgVHSYWgkdFryFDj/L4p0CHVrhoxw2UwUuXSchRIHLnYDQ4n5f78FxWWy6CjZTWXBcBpsui43LYNNlsekq2HRZbHbL4Dh0dE6W9mB16zr5IMG6F+uQ0Y7bk+VAGSkDQxzHw+jmxS9+8WoRmXUgFpG1x4opVQfNFDiMcrSew961jA548jNewlM7GqLApWsHIXLYdAUuwwiRw6YrcBlGiBw2TRfo5Aw2bQ6bLoONzWHTZbDpsuC4DDZdFpsuA47LgEMOHGJNp9uDldd1uilW3Yu1b7SjPVnuuB0HClMYNnI2ehBw1yFw4UBAYmTDGg6Lx4DDXi+mVVo4ZpRTkVHCJoNTp1YZHRa+QYfHz8+2eKfPW7fGputohChwcZ2MEAUusyBEgcs5ENpFp+aw6XLYuBw2XQ6bLoONzWHTZbDpsti4DDZdE+hoIVl7sNxisjteZ99erKnRjqZZGu1omuUgIfZk1aOT89c5zQFcOB+LdRviZFCg4dgc9lix10trOJpOZVxqYOPAyVOrig5w8nNt0KmtsenaQogcNl2Bi8siRA6brsBlGCFy2HQFLsMIkdB5z29YLt4zcCGLTs1h0+XAcTlsuhw2LodNl8Omy2Bjc9h0GWxKGR0Hj9Cp6zrdFKs7ULAb7WhPlqZZGu1ompUhUTo9QqMdNno2fgDhPCz2RoGKYMnptAf2VnEUsqZTeUQjXFwdOBrl8JgAkYVv1qKAc3uks6+AZaobixAFLrMgRA6b0gU6tcDl7AgFLncRGmw/OkpTLPZg6UPd87pON8Xat6BcRzvdorKmWQ4dRhDAw8bNhs4Ig40eNFgQZmGYaZOOu+G/VT7iWNOpOpLJuOSARthUcLTXisfEeg6Pj0Vw4Fyh886fH2goh81UgUvXSQhR4HIjECKHTdcwOrXAZQghh01X4HJ2hAKXa0XIYdPlsOnq0SF3ZLKmWHUvVh3tuAXlqeN26jRLe7McOtqNrgMGmc6wwQOBAAETjWDcwnBGpqLSlbGp4GhaBTiMcnh8LIKDDj/TCp3aJULksOkKXLp2ECKHTVfgMowQOWy6ApY5EVq85zcGKOSw6Qpgrg0hh01X4HJ2hAKXa0DIYUOMdISOm2LlvVgjox23J2tqmqX1HYeOjt/RQYPgw8YufEAAEIQEYIDKFDK6bpegETYOHADksfC4eHwsIjNa42e6+Kusa2y6jkaIAhfXyQhR4DILQhS4nAuhDTo1h40rcLmLUKAxksOm6zB0tBcrT7HcgrKO2alrO1PH7UxNswSPQ6ceOKg1Hn1d8IACOAAFcGRg9G8VFJeQETTCRuDkaZVGOUytWERmtMbPlP8c9KY1Nl1bCJHDpitwcVmEyGHTFbgMI0QOm67AZRghyui8VwBDFh7lsOkKXM6OUOByrQg5bLoCl7MjNI2OTv7M52F1C8r5tIi6J2vfonK3N2u1FpKwUdqVnuFx16OKT0UkQ1L/XenrVLEhTavyKIfHxiIyo7UenVrAMtVJCFHgMgtCFLjMghA5bEyL9/qmQEcFMPdZhAKXa0XIYbPJgUOgo0YWlA8d7bhp1inwuOsQKACEsDg1bkvgaEolcPIoh8fGeg5rU/w8m7/C6rCZKnDpumMQIodNVwDjINpGJxe4XCtCDpuuwOW2IbQXouPRqQvKbrTj1nbyaGdqb5Zb31ktwhpMRgIDUMiLzTkBckj6XgdOHuWAIus5TBH5meofQtwgRA6brsClawchcth0BS7DCJHDpitgmRshD44rcDk7QoHLtSLksOkKXM6OUOAyiJADhzI6bkFZox23tpOnWRrt7NubNQc8QJAXmjNCgkhljLp0XX1/BaeOcpha8dj5mXbQqR2NEAUurpMRosBlFoQocBmGyGFTWrzXNwciymHTFbjcCoQCl7kRmoTIYdMVwDQIOXAoo5M/6mJqtJOP29k3zXLrO+eEhxGH1nsAIQPUIUQZoooM6ft0OwKH+xI4rOWwgMwoh8fP9DH/FY6L1th0bSFEDpuuwMVlESKHTVfgMowQOWy6ApdhhKii894ZnZrDpitwGULIYdMVuNxFaJUDpyZ0RkY7+xaV6/pO3o1eF5YdPN0aj2Kj11oPCFR8HEAZoakyNMKG263gMK1iLYfHy+MH1V10agHLVCchRIHLLAhR4DIMkcNmqsBlFKEVOjULEDlsXIHLtSLksOkKXIYhcti4ApczI+SQqWmKxV6serBgPW4nLyrv25vlFpZH4dEBhDplgg2dKQ0bfV5kdvg4gFSGqCJD+r6MDfeRweEx8Rh5zIza+Bnzn/wZL3DpmgMhsgiRw6YrcBlGiBw2XYFLh9Divb8loMkFLrcCocDFdWUIOWy6ApcTEXLI1PIUS7vP3XE7blHZTbPq+s4UPHmqxV4gNuR98Ezh4wDKCE2l6+r7KzjcL48jj3L4Ofg581/e2PzdMXLYdAUuXTsIkcOmK2BxWYDIYdMVsMyGELXo1AKXIYQcNl2By12EthsAyEFT0zE7mmZptKPjdvIu9DzNmlrfOQWe0VFPxmcKIJUhysAofZ+w4bY1whE4GuUwUuNn2QLHdTRCFLh0nYQQBS6zIESBy7kRWrz3twYkDhtX4OKaAyE6GiAKXFwnIRS4uK4IIYdMTaMdTbPycTt5Udkdu1OnWXVheQqevLis6dYIPB0+DqCMUJeup+8VNtx+BYfHBZA8bn6W+kH4Fp7cFkLksOkKXLp2ECKHTVfgMowQOWy6ApdhhKhFp+bAcQUurpMQCly6rgQhh01X4OI6CaEAxjWIDtXRTl1UrtOsfes7p8AzOurJ+EwBlMsY5X/X9+g2BA73lcHhcfE4edyr9ZwddFyBy1RHI0SBi+tkhChwmQUhClwOQWhx/0AmdzRAFLh03UUo5bDp2sDjgHHl0U7dhZ6nWW5vlqZZh8KTd6e7BWYHT8UnT7scQBmhrgyNsOE2MzjcP4+Fx8bj5LHzs43+7bHt1th0XRlC5LDpClyGESKHTVfgMoXQ4v7fto1O7UoQcth0BS5dRyMUuHTdMIQcMq5906y8N2tqfWcEnnwcD6OeDM/IqCfj40Y/FaF96foZG26f+9EIR+AA5WqUUz4Ef1Xz98emW2PTtYUQOWy6AhaXBYgcNl0By2wI0Q46tcBlqqMQWmPTNTtCDpuuwKXraIQCl64DEHLAuIQO1WlWt75zCDx1d3o+cjlPt+oic4dPnnZVgJQg2peun7HR6Ebg8Hh4jDxmfhaLTu0ohChw6ToJIQpcZkGIApc5EPLo1NbYdN1FKOXAcQUuXTsI0WHoUDfNqus7dWH5GHgY9QiebrpVRz0Vnw6gjFBNuChdn+8XNtwu98H9cf88Fh4Xj5PHzc+UPwT/ooRN1yVC5LDpCly6dhAih01X4DKMEDlsugKXYYTIgEOL+397wEEOm641Nl0HA0RrbLqOQmiNTdfsCDlspgpcXAUgB0yXm2bl9Z18tHJeWJ6Cpx7HoyOXu+nWFD515OMAEkK1Cgzp+oKGNLrJ4PDYeKwXo5wKTtcam66jEaLAxXUyQhS4zIIQBS6HIrR4nwAndzBCCZuuWRFy2HStsek6CqE1Nl1bCJHDpmuNTsnhMlW3vuMWlqfgcbvT8ykTbro1go8b/TiEHEYCRmVouC1hw31xvzwGHg+Pj8e7GuXwOdQH/yVWWmPTtYUQOWy6AheXRYgcNl2ByzBC5LDpClz2IbSDTu0uQuscOF2BS9dJCNHh6JDgYZqVj1bOp0nsgyfvTtcBhGy0gqdOtw7Bx41+MkAVolz+ur4vY6PRTQaHx8jjZgSXPwD/shVC5LDpClimurEIUeAyC0K0g853BC4qYdN17qkYHYxQwqZrVoQcNFMFLl1HIuRgGakuLNc9WhWeurgseNhYgSev8+Tp1qH4OIAyQhmiWr6Ovk/YcLvcB/fHfQscHidoXnwcbMKm6yiEKHDpOgkhClxuI0Lb6OQKNl2zIeSw6UrYdB2M0Bqaqa5kJETnQyePeNwerQpP3aulAwjdOk8d9bi1ng4fB1BGKENUy9cRNMJGo5sKDo8VRPPnUG9K2HRdIkQOm67ApWsHIXLYdAUuwwiRw6YrcBlGiAw2yoPjKth03bEI0RqbrqMRosDFNYGQQ2WkDh42SkYOHOvC+UvAU3en68jlus7jRj11rUf4uJGPAMojICVElECq/67r8/3cjrDhfrjPDA6P24PTtcam62iEKHBxnYwQBS6zIESByzBEGZ33/c7lZRabroJN17nXg+hWIkQOm67AxZXwcaCMJHQED2s34MCGy25nTiXgoyL4gHP+BAx/f0rwcF3gcdOtOurp8KmjH1AQQA4hQTSVrpehETbcF/fP4+Bx8Vj5OS4/DvZdyUEz1Rqbri2EyGHTFY/NZREih01X4DKMEDlsugKXEYS20KlZbLoKNl0Wm66CTde514PoKIQocOk6CSFK+KxzqIwkeNhlDhJsoCzEMsrhrG0+9pM/+8Lfn+JP+vKXNvmrm3WBOU+36qjH4ePWfNz0KyNUMcrlr/M9ug2NbLgf7pP75/EADo87f/70TiuEqEIzVcAy1Y1FiAKXYYgcNlMFMhUhi02XxabLgOOy2HQVbLpuzEiIApepTkTIgTIae7EYtbBhMl3hQDtOnuTT9/iTLhrl8Gd9+RO/qk638qgnr/WM4uMAUkJkKl1X35+x4b64Xx4HjwskmTJuPgo2YdN1FEIUuHSdhBAFLq4rR4gcNl0rdL4rQMkZbLosNi6DTZfFpqtg0zUrQuSw6Vpj07WFEDlstnOgjMSCMRskIwaOceFEST4Ei7+mwB+740/65lGO0OG/gSePetxazwg+GaCKUIaoS9fT93EbGRvuk/vn8fD4eMz5s6d3S9h0XSJEDpuuwKVrByFy2HQFLi6LEDlsuuJ3dk6EdtGpGWxsDpsug02XxaarYNN1MEK0xmaqOUZDgwg5VKZi8Rgg2GDZ68M5S3weDR9qXkc5GRx69KMfvfpa3q2e13r24eMAyiOgjFAuw5LT9/D93I6w4X64Xx4Hj4nHuXi3r44cNlOtsek6GiEKXLpOQogCl1kQosDlWIQWD/ju5eIBDpsuB47LYdNlsOmy2HQVbLpWCJHDpith0zXHaGgHITocHfZMsYGy+MpRvXwsBJ8tPDLK4d9BCZy0dysf0+PwGQFICFWIunQd4nv4/owN98eUj8cCihwSwF+N2PpExtkRIodNV+DStYMQOWy6ApdhhMhh0xW4jCJ0gU7NYeNy2HQ5bLoMNl0Wm6kSNlPNMRo6GiEKXFwJH4eLiwMCQYDpCQfU8TEQ7CJnlMOf82Xx+FGPetRqRONGOYDEIjOLzZwsmo/pmcInLziDTwaoIpQhmkrX5fu4DW6P2+V+uE8eA4+JE1v1J4f5Obfgya0QIodNV8Cyr6MRosDFdTJCFLjMghAFMA4hj04ucJkLoWGIDDZtDpqpCjZdsyNEDpuuBE/KIZNjrxU4MFJgVzMnSWqUw5/jraOcDA7/DUaMcth4Wf/RwYQdPm7NpwJUEcoQ7UvIkEY23L5GNzweHiNACp298NRuy5TsyhEih01XgEOLB3xPoEIOHFfgcptGQmSx6SrYdM0xJTsTQg4bxZHHbKjs9eFAQI1ymHZoF3k3ymF9B5AYDXF9sOI22O2e8WHKNYWPRj95BCSEMkT70vUFDbfH7XNf3C8I8tgAFWSYOu7Cc/FRsBcVbLrmGA1tIUQOm67AxWURIodNV8ByboQ26OQcNl2ByzBEDpuuwOU2TcmOQogCl64TEHLg3HPPPSsAWHTlyN3uQMC6i5w0ymEkxHU5lofvZ2FXp1Do83nyyGcfQBUhQbQvXZfv5Ta4PW6b++A+eQw8JvbIgWoHD4cHbMNTS9h0HYUQBS5dJyFEgcssCFHgcgpCiwc6dGoOm67AZZaREDlsugw2XbOexEoOm67ApWsHIfLoUEWHjZC1j6kDAadGOXydUQ4bLxsze730oWCcrZ7xYeSjDwgDnw6gjFCGaCpdl+8TNNwmt8/9cd88lnsCWY49YkQ2BQ8XD44rYdN1iRA5bLoCl64dhMhh0xW4DCNEDpuuwGUYoWjxwO9dXhS45Cw+ymHTFbjMMhIih81UBZuugxGiNTZdRyNEgYtrACGBw+kObLCMTOqBgGx4+0Y5fE2jHDZgPjyLvUY6dwt89FEZbPD6gDDho9FPHQFlhDJEUwkZ4vu5PWEDeiDISI5Fco4/2gePLqtPZ3x3cth0raGZ6miEKHBxnYwQBS6zIESBS4fQBp1a4HLtCJHDpstBM1XBpmt2hMhh0xW4uCxC37xChw2fvT2cLsBn0OhAQEY5I7vIGeWwi5ypGCMkNmogI45srvho9AM+dfTjEBJEGSOXriNouC1hw/3xGJg+smYleDgkYAqefNHHwl52MEK0xqbrjkCIHDZdgY0Q8uC4ApdZEKLAZRaEIotNV8Gma4UQOWy6ApapTkKIApgGIT4pkPWPeiDgvtMdSAcCMv1iI2V0xEbMbmqdv3VPTGNI+HB/3egHHDJCFaIasOT/1nX5Pm6D2+T2uS/un0VlppAOHn5uB0+97MCTWyFEDpuuNTZdWwiRw6YrcHFZhMhh0xWwzIHQ4oHfF4goh01X4DKMEDlsugKXYYTIYNNlsZkqYdN1FEIUuHSdESE2XA6e40DA173udasNa/R0h3wgIBsqGy6nTQBOPltd3TOBD6OfvP4jiDJGGSXK/6br8n3cBrfFbXM/jLhAienXofC4S/4w/MmOQogCl66TEKLAZRaEKHA5FaHFgzI6NYfNVIHLLAiRw6bLYNP1PuSw6UrYdF0iRA6brsClawchcthsx0YJHM961rOWL3vZy9ab1O5l5EBA9nSxRsIud6HjAh9NuxxAQihDlEHK6d91feL7uS1ul9vnPrk/ru/gYeGcj+1g977g4WdhEX3qUj8Q36JT+ud//uflC1/4wuW//du/LV/xilcsX/va1y7f8IY3LF//+te33XvvvR4ftYMQOWy6ApdhhMhh0xW4DCNEQqdmASIHzVSBy41AiAo2XQcjRGtsuo5GiAIX1yBCTKGe8pSnrDYEXtxTlwwOoxwdCMgoh5ER6yKMHLQ4PZXwuSdGPnndRwDlAMSlr+fvEWTcJref7zPDw54t4GEayPlZGR4OiBQ8+y678Gz3xCc+cfmc5zxn+S//8i/Lf/3Xf12+/OUvX7761a9eYfLmN795+aY3vcnG18hd+He+3wJEJyNEgcssCFHgMgXR4kHfH9DkApdhhMhh0xW4DCNEDpsuB81UBZuuFULksOlaY9O1hRA5bLoCF5dBiD1MAMJGwbvu6IXv0YGAjHJ0ICC72TkYL2/oU2V4VB4B1fhaLn9fjtvktt19cjvAw5QSeHi8HTxs2ACw71Kh4ffzpCc9aTV6FDaMIl/1qletpq9vfOMb19+5uQiZDpruwvX3TsksQuSw6QpchhEih01XQJMR2kXHFbjMghAFLrMgFFlsuhI0U904hCiQKbHuwprNU5/61OWLX/zi1RD/kEs+3SEfCOg29KnAQfhUgGrCpFZvc9/Z9cDDuk8HD4cMvPKVr1y+5jWvGYfnfo9Yjf4e//jHL//pn/5p+dznPnf5ohe9aPmSl7xkBfoht3XsZXtN6IoR2guRw6Zp8aAfCEjIYdMVuAwjRA6brsBlGCEy2HRZbKZK2HQdhRAFLlOdiBCnNzz2sY9dbRy8Ax9zyQcC8gHobgMfrUIihEahoXxy60X92fVMyQQPBxIKHuDksAFGJsDDmssoFoxunv70p6/WbSo4TF3nBCdfdvHJBS6zIESByzkQ2qBTc9hMFbjMghA5bLoMNl3vSw6broRN1yVC5LDpWmPTtYUQOWwuYm8NYPCOzILmsRtDPhCQE0Tdxn1oDhjlrk+72GznvocED8f0AA/naz3hCU9YPuMZz1i+4AUvOAqe6wYnXzw6tcDFZREih01X4DKMEAmdBwcwOQsQOWimClyGESKHTZfDZqqCTdfBCNEam6nmGA1NIMT6CxvX8573vNVGceyF9QQOBGQdxG3UpzYFDdWjqy/y8HQL3MDDbnbgAWL2yPG7eeYzn7laj6nw7Ftz4XuAXAvGgOPWcK7qcrbd9RYgcthMFbjsQ2jx4B8MbFTgMowQOWy6ApZbiRA5bLoSNl1zjIbW+HDwH2sxvJvzbnzshY2Pd28WXDmSmYVpt1Gfs81R1BWaqabh4d+YGrK3iuNyWBhnbYa9ThkepqCj8ORF40PXyua4bI4VyjlsugKXYYTIYdMVyFSEttFxBS6zIESByywIRRabrgTNvuYYDR2NEG3jw3SItQemD2xEx14Ah3dwDii8XKAue8dOhSgfNT2dw8YHMq4peDicIMMDJPvg0cLxvutd5eXyeKGTEaLAZRaEosWDf2iNy2jxIhxGiBw2XYHLMEJksGlz2ExVsOm6QQixS5u9TvsOBNx3YUMCHDY+PihrM3ULaHIGoe12jxvqC1yGESKPDjl0SJ+SyPljHImdd30DD1OmQ+C5zmlVd8kHKt5YhC7QqTlspgpcZkGIHDZdDpuJLDZdBZuuOaZkWwiRR4eN6clPfvLQgYBTF0Y5bHSsCeW1ot0CltkQIodN1+HwMA3N8LBAfCg8N2mko8vmmKLAxWURIodNV+AyjBBldN4vkFEWIHLQTBW4DCNEDpuuwGUuiCw2UyVsuo5CiAKXLoMQay9MGZ797GcfdCBgvWiUA1osunpspgpcbghCDh2O0+HgQBbIOVpb8IA1v7sMj9Zsbjc8tcBlFoQocBlBaPF+PxzgUOBSswAph01XwDIbQuSw6TLYdD2AHDZdCZuuS4TIYdMVuDRx3AnnRnEgIBsNG8qxF41yeOff3juWc9h0BS7DCJHDpsths90+eDgOiY/t4LyzKXj2XW4XPLnAZRghcth0BTAOoQ06rsBlFoQocBmGyEEzlcOmq0Az1cEI0RqbrqMRogt0OCfqMY95zOpAQHb/HnvRKIddwXWB2uOjHDZTBS6zIERj8AA1a2Ac+Fjh4bQRjuLWbvGbspfq0MviPR6xzmEzVeAyC0K0F51a4DKMEDlsugKXYYTIYdPlsJmqYNN1BoTYta0zk3XAmXbd5pju5P9PwMNfdOAgPg4E5BygUxY3Ncrhtup60WEIkcOmK3AZRogcNl374eG4Js4v41yzf/zHf1ydIOvguYkLx/suG3hyDpqpApdzIrR4vx8JQHIOm6kCl2GIHDZTBS6zIBRZbLoKNl0rhMhhs4kXdf4IBF7Y+dydLkYi3TCef2fjOPaiUQ57vXbB6QpcZkGIApczj4YcPJyjxicpOngYOQI5bwS3Gp77BTTKIkQOm67AZRghqug8JKBROwAph01X4DKMEDlsugKXYYTIYNNlsZkqYdOVEMrvoGCTPwaBUQuo5IuAUaOXQ66rC9/D/TPKYU9OXqA+rMBlGCIHzVSByxlGQx08fE4QH3XB2hjwMF1ljexOgGdxv69fF7jULEDksJkqcBlFaPGQH92gU7MAkcOmK2CZDSFy2HQZbLoeSA6broRNij0jHITGiYLsgmbI/tKXvvSgQ++PvYzeLuCwIfHYVutEW3vGag6brsBlGCFy2HQFLEcgxJTKwcNBkBkePqI1w6MP5dKbxO2FJxe4DCNEDpuuwKVD6AIdV+AyC0IUuNzhCLErlhctC5Mc9VrBuSlHtHL/bECMcjgrPS9Qb2UBIofNVIHLLAiRw2Y7RjUcndzBwycL8jnQfCwI8PA7edrTnraaDt+Z8NQCl1kQogCHPDiuwGUYIXLYdAUuwwiRw6bLYTNVwaarQYiDzvRi7cCZe4RzyEXTquc///mBS907pgKXYYTIYdMVuAwjRA6brm1wBAsHCPL5QBUdYmE+w8MR3hmeenb5bYJncb9viBw0UwUuwwiRw6a0eMiPBSDKYTNV4DILQhS4zIJQZLHpKth0BUJaiBQ42jNVPwpBC8bXfdEoh8fDNNCD0xW4zIIQBS4zIAQqQofDDPiLGBUd4qNBONiSr/P33QUPe/V4A8nP5W27XMBTc9h0BSwHQeTQef/AhrbwyTlsugKXYYTIYdMVuAwjRAabLovNVLvoAI6O9WAxlheoPjv3Jn0UQr5olMPj3D5eiBw0UwUuNwIh8uBwUCDocGAgJ3+CCsfoVHSowsOH2z/ucY/bgoff300YrR5yyZ+66AEih81UgcshCF2iU7MAkcNmqsBlFoSoQjOVwabrQeSw2Y0PuOKFywmEnIbAZ7Ww+MgRvUxZ6rSKjfymXMCPd+t67NAuQMph0xW4DCNEDpuugOVAhICERWSmTxybk9Hhg+crOsTnCPH8cuY+x0JleHge+f3dPnT0USCBi8siRA6brsBlCqHF+/94IFMLXO54hKhg09UgxNyfNRw+HIoXpBaO2VPFHg83yrkJ06p84bFuH0cUuMyCEAUusyBEDpuLAIRRDovIfD4yx+UIHf1ZZUapFR2q8PCB9VqTu43ocNnAUwtcZkGIAhvl0akFLsMIkcOmK3AZRogcNl0OmqkKNl2BEEN0Db15wbImwiiHBUdGOQzBWcupoxxepDfthbqNTi1wGUaIHDZdgcswQuSw6dpFh1FORqf+LXdAqegQ38PCc11Ezs+ne04Xi8VO7jJ6vXNdPDiuwGUYIXLYmBbv/xOBCDlspgpcZkGIApdZEIosNl0GnIgXoc5S5gOhNMphLYeRA2skfIiW/jRJHuVMvUiv8+LB6QpcZkGIApczIwQcGunwyYqgw0Iy6DBS5Y0DdHgeWZPjLz7wPTyXnADKG0g+epznU+djTT2fo5iMXu9cl8V7flPKYdMVsJwDocVDhY7LYdMVuAwjRA6brsBlGCEy2HRZbPpYiNRwW6Mc1nKY63PEcT6XSi/SvMdq6kV63ZfR0zh2C1yGESKHTVfgMowQ7YKjNL1iTQd0OBCQTxLUn1dmLxUjVp5L7Qxg1MrzyVRZbyL5sIf8nNbLKCaj1zvnZfW5Q1v4KIfNVIHLoQit0KlZgMhhM1XgMgtC5LDpMthMZbBROjNZazkck8O7o5ta5eF4Xni8qehw2cBTc9h0BSw3BKGMjqZYTI2ZLnG+FX+5lCOQGbHyXPIGwu5xjiLn+eQ4K9bnGLmyU0DPad4T2T2Xo5iMXu+cl/yBZ5dZhMhh0xW47ENo8dCfDGhygcttRWjm0RBDc0Y5rOWwi5x3R4bknFfFi5R3Robj3V4rXqC3YTfr1qkcFiBy2EwVuMyCEHlwqKKj43R4Hnnz4A8J8lyyB4tRKx9cxvPJmwgj126KpU8X1IKyu4xiMnq9c17qx3/MhxAFNBmhXXRqgcswQuSw6QpchhEih01XwHJGhLSrlY+65IXKkFyLjzoYUOs5FR3N/28POjrCOnCpWYCUw6YrcBlGiBw2XbvoAI7Q4TQIjsFhXY5pMrvNtZisdR2mWDyndYql57UuKPO81ssoJqPXO+dlF51a4DKMEDlsmhYP/amAxGEzVeAyC0IUuAxD5LCZymHTtQGHFytrAewi19SKFylHH2s9J6OjoTgvTn36XJ7/34ZLPrVju8BlFoQocDkzQhUdDvpjxMrJnbyBaDFZ6zp1isUxV6zT6c2EaXMd7XTTLIfJaHNfeD16bLoClnMhdIGOy2HTFbgMI0QOm67AZRghcth0OWx207sj6wDazcpwXC9QLSJXdPKfNblt6HDx6NQCl2GEyGHTFbgMI0QeHKHDwYE8l9ptrsVkreuwRscbCaNXpljsxdKbSTfa0Z4sN81ymIw254XXIFAu3utbSg6bqQKXYxBaPCyAyVmAyGEzVeAyDJHDZqrAZRaEogIOL1gtPup8HYbjWkTW3g7tuXIjnX1z/5t6uTiplRw2XQHLbAhR4DIIUUZHAQ8H/PF85nWdOsXieWW0w14sLSjX0Y7eVJhm5V3oujhMRpvzwmPk8ebPGdoFSDlsugKXEYQWD/vpwEYFLq4rR4gcNl2ByzBCZLBp4oXKkFyHzrMGwMIj6GgXK++IbqSjNZ3bOtLhsoEn57CZKnC5BoQqOIo3EUauTJeZYvFGoimW9mK5BeU62sl7supzTA6T0ea8bNDRGfiBS80CRA6bqQKZitA2OrXAxWURIodNV8BywxHiBap1AIbjWkQGHV6c+biOjE5dSNZxOm7ef9MvWye2WoCUw6YrcJkNogBnjVDFRs9nPjI5T7HyXqy8oKzRTl7b0ac/6ridfMAg8NAoJqPXO9eF1yBI5o/92C5wGUaIHDZde9FxBS6zIESByywIkcOmaxcdFh7Z2wE6Oq4jo6O5v3aZu4+z0N6r23bZnG8WuNQsQOSwmSpwOTNCFRw9n6zr8JzmKRYjWNbqOOBTC8p5tJP3ZOm4nbyorDcXLSrTKCaj1zvXBXR4rB4cV8ByToQWD/uZgKPmsOkKXIYRIodNV+AyjBA5bLocNpv0AtUaQEZH06u8pqMXZXfY/NSu1Zt+ySe5bhe4DCNEDpuuwGUYIRpHR/BoiqXzsDg6mTeUPNrR2k7ek6XnOE+z6rE7PNejmIxe71wXXoO8Jhfv/W0ph81UgcuxCC0eHsjkTkaIApdZEKLAZRaEIoMOL052sfLidOho92o9bF7oaLFRQ+/bOMXi4tGpBSyzIUSBywEIVXD0nFI+D4u9knW0U9d2tCerLiprmpX3ZgHPas1kEJPR653rAjqMzC4/b2gLn5zDpitwGYVo8fCfDWwocHFdOULksOkKXIYRIoONSS9O5v8ZHS0k827I7lVelJzsWQ+b17k6dTFZ6Dh4Rl98o9c75+Xi4z3IYTNV4HJNCE2hoxEsRydrQVmjHT2/7MliCs0bC8+xFpXzNEvw8AZT4XHP0ynP+ykXveaINz8e78VZ+IFLzQJEDpupAhiH0AadWuDisgiRw6YrYLnBCOkFquM68p4OHcGqIXie97sDyfK6Tj6I7JQX3+j1znlZnWF/CU/NYdMVuAwjRA6brsAmIZSxUYCj/9UbihvtaE+WjtvRorKmWVq/Y1Sb13fywrJ7nhhlZABWv1tzvXNedF/cN+AAIq9PffTHdoGLyyJEDpuuNT4enK7AZRaEKHCZBSGq0Ey1QSe/MHUEa/7gp/pxCNqD1Z2roylWfsHly+iLb/R657zUj/c4H0IUuMyAkJ7HLt5Q6mgn78niOa6LyjzPeW9WBw9vNO55YsNXV40Ob3i88fFa9OB0BS7nRGjx8J8LPHIOm67AZRghcth0BS7DCJHDpsthsym/MDUM591QC476DJbRj0Oou87zCy5fRl98o9c754XHm8+2Pwwhcth0BS7DCJFHh/JzWWNdR7vPOQ5Lox0dt5MXlfM0K6/vaGSrheUKDyMennfecDS9pu41MMdFyHH/PBZej4v7f0fKQTNV4HIKQosPCGjUDkDksJkqcJkFIQpcZkEoMuho/q+DA/NnsHQfh8DQW+s6dYrFCzAfy1Evo5iMXu/cl210agHLbAhR4HIEQhmZGs+tdp/nPVl85AU7DPI0q+7NyrvRtbDcwaO9WtqDeZXocB/cF/erUQ6Pb+tzh7YAUg6brsBlGKFo8QE/v0GnZhEih01X4DKMEDlsugKXYYTIYFOqL0qho2M69BkseTE5n6vjplj1w5/04qsvulFMRq937ovHZqrA5ZoRys+ni9GOplm8sXD2uRaVeZ51egSj2ry+I3i0sJzhyVMtRhWaXtcRz1XAI3R4zXH/vA5X6zmrs/ADF9fJCFHg0iF0gU4tYHFZgMhh0xWwzIYQOWy69qOT93JMfRxCnmLlvVg6JaKOdvTCyxeHyWhXcVk8+IdNDpuuwGUYIXLYdAUyBqH8fLp4jjXNqovKTLPysTt5fUfwMLqt8OQ1Hk2vdQCh9myBAK8B3nzy6EcQqVMuAofb5zXH/fN4mBLqoz+2C1xcFiFy2HSt8SGPjitwmQUhClzmROiAKVn3gsx7sLSuk8/V0Xw/H0Smd7462uEdR6MdXhS6OExGm/vCCzh/3IcHiBw2UwUusyBEY/BoL2VeVGYqzaiWN5i8vsObTIZHe7QyPHlxWc89bzp5usVrAHzAQABVhDJEKoO0L64vcLhPHgOvRw9OV+ByboQWH/ALAYhy2HQFLsMIkcOmK3AZRogcNl2BywRC9QWp0c6+c3U07NYioxaU62hHRyjr3S6Pdhwmo8194UW8Ofs+cKlZgJTDpitwGUaIHDbb1ee0pudYi8pMpRnVsr7Dc53XdzgoVMfvOHjyXi2t6+nI5Tzd4s0HCPJCs3IIVYy68nW5HV5n3Bf3y+O4GOV8V8pBM1XgMowQOXQ+MLChLXyUw2aqwGUWhChwGYbIYTNVjw7xgnTn6uT5vvZu5AVlvfjyaIcnX8NsXmC8OPTO5DAZbe4Ljw8888d+bBe4zIIQBS4nIOSe05qbZvFc1/WdDp68uCx49KaT13nyqMfho9FPBkgJolq+jr6P2+E2uX3uj9cfI28e49bnDm0BpBw2XQHLoQhdolOzCJHDpitwGUaIHDZdgcswQuSw2c29ILWu002xeBG6D3/ScFujnbxXQ0NsXii8QzEUHsVk9HrnvPD4eNweHFfgciMQosPg0d4sre+wjncoPDqAMK/zaNRT8WHKlfERQEJIZYxy+rq+R9iAGvfB/XH/YMio/OKjPwIX18kIUeAyhdDiA38xkMkFLC4LEDlsugKWgyBy2EwVuJyIkHsxag9HnmJpt2peUOYFWE8QzEeuapitaRYvMl4svFsBj8OEEUa9uOvNfdmg86Mph01X4DKMEFVopgpY9kDknteuDA8Hhe6DhzUe3nB47nUcDzsUtM6Tp1sa9TDyyPho2qU1nwxQRqimr3N9oBE23B5vcNwPrz0ew8UoZ/N5Q9sFLi6LEDlsugKa3C46rsBlFoQocLkRCJFHR3P+PMXiRVg//KkeMp/3ZPFuxztdnmbVvVkOk5uEDo959bEfW/DkHDZTBS6zIESBTUHIPbcujXi0sCx4mFZXePLisuBhtKt1Hk23NOqZwgckNPoRQIJEEOX4d65DXJ/vEza8ufEmBzjcL4+Dx+jB6Qpc5kBo8YG/FIiQw6YrcHFZhMhh0xW4DCNEDpuugGUPQt2LkLQXixehFhm1oJxHO1MnCNZpFi8c3rGAx2Fyk9DhhazPGtrKAqQcNl2By2wQHTba4U3G7dHKIx4Wl/PudN50OHyC51/TrTzqyWs9GR9NuwCI33EGCESI14owUvo3vs51hQ23w21y+wIHDBfv+z0lB81UgcswQmTAocUHBTh0iY9y2EwVuMyCEAUusyBE4+hoisVBZG60w3A778nSu12dZnXwOEyuCx3uVwEOa0/5oz+2C1yGIXLYTBW4nAkh99xOBTxuxCN4tFcLeHjT0ZHLWufRdIvXQV7r0ZQr48PrQqMfAcRrhASRi69xHa7P9wkbbpf7ATxei7w+6+cN7SJEDpuugOUYhC7Rqe0gRA6brsBlGCFy2HQFLsMIkcPG5158gMP/aopVRzt5bScfQOZOENSxO3V9B3gcJmzwGQBy1zv3RffF/TP1YxrowXEFLsMIkcOmK3AZRoiOh4fnfQoePfe86TDa5Xy8vM6TRz0OnzztcgBpFAREitcN6b/5Otfje3htcTvcJrfPGx73zePJH/mxg486GSEKXPYhtPigXw5kVMDSdTJCFLjMghAFLmdAyL34lKZYbrSTj1x10yxeaNqb1cHjMGGXOhu+ukp0uD8WuJn6MSJbPOTHUw6broDlBiHkntuuDh7tTtcBhDz/eZ1Hox6t9Th8NO3S6KcCRLxWgMjF17gO1+f7hA23y30w0uIx8Bi30akFLi6LEDlsugKZ2jY6rsDFdUchRPvRobpnI+/JYqitT53TojJPvPZmaX2HFwgvGs3l64iHqRZ7tZjSaM+WdquDwdyXDA73z+PhMV5+5tAWPsphM1XgMgtCFLhMIOSe16kyPHl3ug4gzAvMvAbyqEdrPbwJOXwYCfPayAAxAhJCjI4VrxvSf/N1rsf1+V5ug9vjthlpAw4QLh7wfZHDZqrAZRaEosUH/Uog4rDpClxcFiFy2HQFLsMIkcOmK2AZRMi98BQvPkY7vOvpyNV8AJkWlTXNqus7vDh4sfDuJHi0gAg8WuPJe7UY7VwHOtw3+LF3hKH8JTo1ixA5bLoCl2GEyGHTFdgUhNxzO5XgYV1PBxAy2s0LzEy38qhHaz0dPtrTxZtSBYgYHQuiGv/O17keoyZhw20yrWO0zf1fnnG/gqfmsOkKXFwWITLYqMUHg07NYTNV4DILQhS4zIIQHY4OLz7e8fKRq9qdWj+HRdMsre/wwuJFwoumwqMRjxaXtTtdIx7BcxXoaJTD/fIYwJChfP74j02ByywIUeAyC0J0ODyU4dEpE5pu5VGP1no05cr4aNoFDrw2MkDAIYQEEQFLjn/j61yX7+P7wYzbZaTNiCt/xMdO14nQCp3aDkLksOkKXIYRIodNV+AyjBA5bLr2o0N64THUzrtT8zQr783S+k5eWK7w5DUe4GE3KPBwPEaeboFPHvlo9KNOvWiUw31olMNIzIPjClyGESKHTVfgMowQOWy2c8/vvvKohxEv0y19LAavA6318CbE6Ff48GbEtEtrPrwpMfoRQKCRESJGyTV9jevxPby+eJ1xm7zuuE8eVz3j/qLApeskhChwGUFo8cG/GtBQwNJ1MkIUuMyCEAUuZ0TIvdByGu1oF3qeZuW9WRke7cno4MmLy+wCZXco8GidJ+PDtEvrPRkhJYwcSvvK4HC/AHgxyvnJlMOmK2C54Qi553gkXgeCh1GvRj3CR1Mu8OE1wbRLaz559COAeJ0IIQIT4rWj/098jetxfUZOIMbtcduMshh9e3BcgUvXDkLksOkKYBxCG3RcgYvryhEih01X4DKMEB0PD6OdfIKgjuHg3QZ4tLA8BU/eq6XjeOo6DyMO4aNpF/hkgJQgqhh16bp8L7elEQ7gMPVj/eDyM4e28FEVmn0FLrMgRIHLgQi553e0/AakUU+ecgkfjXzAJ49+BBCvk4wQAUpO/851uC7fw/dzW6wlcT/c9+pE1wfmHDZTBS6ukxGivejUAheXRYgcNl0Byw1CyL3Aanqny3uz8voO72zak+HgYU5e92rpWA2t82jU0+GTRz8ZIZUxyuXrEN/H7XCb3Af3x/0zCmPIf4lOzSJEAmakwGUYIXLYdAUuAwi553c0waMdDIx+Mz552sVrQ6MfoAAgXicZIaZhBEY5/o2vcz2uz/dyO9wmrzvuM59dv9UWQuSw6QpcXBYhctikFh/8awFHzmEzVeAyC0IUuMyCEDlsNrkXl4u5PS84Xmx5fedQeNgTkQ+T13Qrj3q01tPhI4BUhiiXr6Pv43a4PW6b+wEcHsPFKGfz0R/bBS6zIESByywI0Tzw6LXAG5EA0rSLkTCvDY1+BBAjoIwQaTQkkEj/zde4Lt/DyIbbAjbuN59nZuHJXSdCiw8JaHI7CJHDpitwGUaIHDZdgcswQuSw6QpsCkTuxeXSiIcXGS+wvLDcwcPcnBEEC4PsgWCPhI5UZZ0nj3ocPhr55DUfAZQRcuk6ggbAuB1AY1QFONwvoy9g9OC4ApdhhKhCM1XgMowQOWy6TseH14ACnwyQpuC8KYGPRj8AxAiFvV4EREACRi6+xvV4bfH93NbF6CZjM1Xg0nUSQhS4jCK0+JBfD2wocHGdjBAFLrMgRIHLMEQOmz734uoSPHlhmRcI8OjYDcGT92oJHnaVss6j6VYd9bDWw5Qr46Npl0Y/wJFHQC4hI2j4Pr4fyJjKcfuAA3w8lsVDf7rksOkKWGZDiAKXWRD6Mfscj5TxUQDEyAd8eH0IIPZ4AYcg4s2K102Nf+frXI/v4Xu5HV5vq9M7HpRz2HTF93ftIEQOm67ApUNog04tcHFdOULksOkKXIYRIg+Oci+sLl5cgke7UAWPRjxaXAYe9liwJ4LdnuwC1XQrj3oYbVR8NPLRmo8A0ghICAminP6d63Bdvo/vBzJQ4z64Pxa5QTF/1tAuQOSwmSpwuSUIuef42MCnrvsoEGLEov9V/DfpenwfcPEa4/byOWVbHY0QBS6ukxGiNT4eHFfg4rIIkcOmK2A5CCKHTVfAcgJC7gXUpRcWLwpeLBkehsfAw7wceFgYZFcneyXY/al1njrq6fDRtEujHwACDiEkiFx8jetxfb6X2+D2uG3uh/vlsWzOvg9cXBYhcth0BS7DCJHDpitwOQNE7rk+tvwa0eiH/83xNd7ANErKcRv5yOqL1th0bSFEDpuuwMVlESKHTWnxIY8MPJTDZqrAZRaEKHCZBSEKXA5AqL5wptILRSMeTbW0uJzP0WFvBLs+tc5TRz2s9XT4aM1HAGkEJIQEUU7/znW4bj0zmWkd98f987jyR37sFrjMghAFLrMgRIHLkQi55/uYKiYgRPx/vqbq9+1i0xWwTHUSQhS4nILQ4kMDG7UFkHLYdAUuLosQOWy6ApdhhMhh0xW4TCBUn/x9OXiYl+/7SASNerTW0+GjNR+gEEDAIYQEkYuv6cxkvpfbADOmU9wP9wuCi4f9bOSw6QpchhEih01X4DIbRIHLMEJ0Xnxq7nr5QMbFg5XDZqrApeuqEVp86G9s0HGdhBAFLrMgRIHLLAjR+eDR7nQWA4FHh8lrgZnpVh71aK2HKVfFR2s+IJEBYgQkhBQYKf0b1+G69cxkbp9RFvAB4gU6tQrNVAHLbAhR4DILQuSw8bnn/hxd7ElL2HRdIkQOm67ApWsHIXLYdAUuUwhdoFMLXLquFCFy2HQFLsMIkcOm6zh4NH8HHhYH2QPRnRzIqEdrPQ6fvOaTAWIElBFi5OLia1yH6+qkQZ2ZzOiK++Rx5M8XuuxkhChwGYbIQTNV4DKMEDlsujw4Lvc6GC3vtu9bY9N1NEIUuLhORojW+JBHpxa4dO0gRA6broBlNoQocDkTQu6Fsq8Mj47XYIFZ0y2NerTInM/NqScGAoNOCgQKAGKEkhECEgKVGv/Odbgu38dtcHvcNvfF/dcz7i1A6mSIApdhhMhh0xWw3ACELtpeG2p7P+Ww6Vpj07WFEDlsugIXl0WIHDamxYfGUHqVw2aqwMV1MkIUuNxQhBwsIwkf9k6wzqPplkY9WuvRlMudGMjURwvOAMQISAgxChJELr5GXI/v4fu5LW6XqR33uwtOV+DiOhkhClxmQYgCl1kQiiw2XQmaqW4cQhS4nIrQ4sMCHHUJkKrQTBW4uCxC5LDpClyGESKHTVfgMowQHQ8Psc4DPJpuadRTp1wgwMhD5+QAAwvOjEgEUD0zGUiEUY5/4+tcj+vzvdwOoHH7jLQWD//5ksOmK3BxWYTIYdMVuAwjRA6brsBlGCEy2HRZbKZK2HQdhRAFLl1zIEQWIaro1E5CiAKXWRCiwGUWhMhhs50DZTQ36slTLuGjc3JY89HoRwBpBCSEmIYRoNT4d67D9YGLERTYcNvcF4+hnmV/GkIUuMyCEAUusyBEDpsug03XQ8hh05Ww6bpEiBw2XfF8d+0gRA6broDFleFZfFjM47cKXLp2ECKHTVfgMowQOWy6ApdhhMhh03V+eLRni3UeHaGa8cnTLkYijH4EkEZAQghEBFGNfyeux/fwvdwOt8vxQzyGzflmgUvXHYMQOWy6HDZTFWy6DkaI1th0HY0QBS5dJyFEa3jULjq1gKXrRiNEgcssCEUnwsNoh8BHC80a+WjapRMDNfoBCo2AhBDTMEFEjIhy/Btf57o6M5nb4/bzya2+wKVrByFy2HQFLsMIEbiMFrDcSoTIYdO1xqZrCyFy2HQFLl07CJHDpmnxYb8dgDhspgpcXCcjRIHLjUCIDDYmh8pIGR6laZfWfPLHIgggjm4GIXdWMqMi0n/zda6rM5PB5mI69YvrHDZTBS6ukxGiwGUWhChwmQWhyGLTlaCZ6iiEKHCZ6miEKF43rkMRWnx4oEMrfHIOm67AxWURIodNV+AyjBA5bLoCljMi5GDZl+CpZXw0+hFA3VnJTMn0//ka1yGg4Xu5ndV0qpxftgFIOWy6AheXRYgcNl2ByzBC5LDpClyGESKDTZfFZqqEzVS3GiGq6NROQogCl1kQosBlFoTIYdN1Hni6HD6ks4+BqKav6Xp8H4ENt3dxhHXAMtVJCFHgMgtCFLjMghA5bLoMNl3vTw6broJN122ZklV8Fh/+O4GMCly6dhAih01X4DKMEDlsugKXYYTIYdMVuBwBkUPkmICirvuw+FzPTubreZqm+H5uJ5/S4Qtcuu4YhMhh0+WwmapgM5XFZqqETdccU7KTEKJApraNTi1g6ToZIQpcZkGIApdhiBw2UzlsfBWRYxM+CmQEDV9T9fsWH/jLKYdNV+DStYMQOWy6ApdhhMhh0xWwzIZQZLHpMth0zTEaukSIHDZdgUvXDkLksOnai44rcHHdUQiRw6bLg5OrGBxbBqaDhuqR1JcdjRAFLq6TEaLAZRaEKHCZC6EbMyWjNTZdRyNEgYvrGIQWH/G7gUbNYdMVuLgsQuSw6QpcZoMoYLkGhBwQ5yofQ7RpjU3XFkLksOkKXFwWIXLYdAUuwwiRw6YrcJkNIofNVAWbrluNEFV0aichRIHLLAhR4DILQhS4zIJQZAAih8chLT7oVyOHzVSBS9dJCFHgMgtCFLjMghAFLnONhiw2XQWbrhVC5LDpWmPTtYUQOWy6AhaXA2jxEb8X0FDg0rWDEDlsugIXl0WIHDZdgcswQuSw6QpchhEig02XAWi67b1kthVA5KCZKnDpulKEyGHTFbjMBlHgMowQGWy6LDZTJWy6jkKIApeukxCiAMYhtEGnFrB0nYwQBS6zIESByywIkcOmy2DT9cHksOlK2HRdIkQOm67ApWsHIXLYdAUstxIhcth0GWy6HkoOm66ETdclQuSw6QpcunYQIodN1xoeD05X4NJ1xyBEDpsuh81UBZuugxGiNTZdRyNEgYvrZIQocBmGyGEzVeByIxCigk3XwQjRGpuuoxGiwMV1LEKLj/z9QCPnsOkKXLp2ECKHTVfgMowQOWy6ApbZEIosNl0Fm64VQuSw6Vpj07WFEDlsugKXrpMhClyGESKHTVfgMowQOWy6HDRTFWy6ZkeIHDZdgYvLIkQOndrRCFHg4joZIQpcZkGIApc7HiEKXLpOQogCF9fJCFHgMgtCFLjMglBksekq2HStECKHTVfAMtVVIrT4yD/YRae2hRA5bLoCF5dFiBw2XYHLMELksOkKXIYRIoNNl8VmqoRN1xwI0ZUhRA6brsBlGCFy2HQFLsMIkcGmy2IzVcKm6yiEKHDpOgkhClxcG3RcgUvXSQhR4DILQhS4zIIQOWy6DDZdH0IOm66ETdclQuSw6Vpj07WFEDlsugIWlwWIHDZTBS6zIEQOmy6DTdfDyGHTlbDpukSIHDZdgUvXDkLksOlK8HhwXIFL15UiRA6brsBlGCFy2HQ5bKYq2HQdjBCtsZlqjtHQSQhR4HLHI0QFm66DEaI1Nl1HI0SBS9exCC0+8g8DDuWw6QpcunYQIodNV8ByKxGKLDZdBZuuFULksOlK2HTNMRraQYgcNl2ByzBC5LDpClyGESKHTZeDZqqCTdfsCJHDpitw6dpBiBw6HxXY5I5GiAIX18kIUeAyC0IUuNwmhGiO0dDRCFHg4joZIQpcZkGIApdZEIosNl0Fm64VQuSw6QpY9nU0QhS4uDqEdtCpbSFEDpuuwMVlESKHTVfgMowQOWy6ApdhhMhg0+awmapg0zU7QuSw6QpcXBYhcth0BS7DCJHDpitwGUaIDDZdFpupEjZT3ZYpmeBZfNQfBS61wKXrJIQocJkFIQpcZkGIHDZdDpuJLDZdBZuuOaZkNxYhClxmQYgcNl0Gm66Hk8Omq2DTNcdoaAshcth0JXjIo1MLXLquFCFy2HQFLsMIkcOmK3CZCyKLzVQJm66jEKLApeskhChwueMRooLNVBabqRI2XUchRIFL19EI/fTy/wepOKVWiJemWgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 2,
     "metadata": {
      "image/png": {
       "width": 150
      }
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import Image\n",
    "Image(filename=\"assets/ch4_mol.png\", width=150)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "分子对应的 [输入卡](CH4.gjf)、[输出文件](CH4.out) 与 [fchk 文件](CH4.fchk) 在链接中可供下载。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                     1                      2                      3\n",
      "                     A                      A                      A\n",
      " Frequencies --  1675.7695              1675.7695              1675.7695\n",
      " Red. masses --     1.1714                 1.1714                 1.1714\n",
      " Frc consts  --     1.9382                 1.9382                 1.9382\n",
      " IR Inten    --     6.6431                 6.6431                 6.6431\n",
      " Raman Activ --     3.8439                 3.8439                 3.8439\n",
      " Depolar (P) --     0.7500                 0.7500                 0.7500\n",
      " Depolar (U) --     0.8571                 0.8571                 0.8571\n",
      "  Atom  AN      X      Y      Z        X      Y      Z        X      Y      Z\n",
      "     1   6     0.12   0.00   0.00     0.00   0.02   0.12     0.00   0.12  -0.02\n",
      "     2   1    -0.60   0.00   0.00     0.00  -0.08   0.11     0.00  -0.60  -0.01\n",
      "     3   1    -0.60   0.00   0.00     0.00   0.23  -0.49     0.00   0.01   0.29\n",
      "     4   1    -0.12   0.28  -0.20    -0.16  -0.17  -0.53     0.30  -0.43  -0.05\n",
      "     5   1    -0.12  -0.28   0.20     0.16  -0.17  -0.53    -0.30  -0.43  -0.05\n"
     ]
    }
   ],
   "source": [
    "with open(\"CH4.out\", \"r\") as f:\n",
    "    while \"and normal coordinates\" not in f.readline(): continue\n",
    "    for _ in range(15): print(f.readline()[:-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们这份文档的目标是重复其中的红外光谱强度 `IR Inten` 一项，并依据这一项对红外光谱进行绘制。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 简正坐标的导出"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这里的代码将会给出归一化后的简正坐标 `q_normed`。这是 [上一篇频率分析文档](freq_1.ipynb) 的主要内容，这里简单作一个回顾。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from formchk_interface import FormchkInterface\n",
    "import numpy as np\n",
    "from functools import partial\n",
    "import scipy\n",
    "\n",
    "np.set_printoptions(5, linewidth=150, suppress=True)\n",
    "np.einsum = partial(np.einsum, optimize=[\"greedy\", 1024 ** 3 * 2 / 8])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# https://docs.scipy.org/doc/scipy/reference/constants.html\n",
    "from scipy.constants import physical_constants\n",
    "\n",
    "E_h = physical_constants[\"Hartree energy\"][0]\n",
    "a_0 = physical_constants[\"Bohr radius\"][0]\n",
    "N_A = physical_constants[\"Avogadro constant\"][0]\n",
    "c_0 = physical_constants[\"speed of light in vacuum\"][0]\n",
    "e_c = physical_constants[\"elementary charge\"][0]\n",
    "e_0 = physical_constants[\"electric constant\"][0]\n",
    "mu_0 = physical_constants[\"mag. constant\"][0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "class FreqAnal:\n",
    "    \n",
    "    def __init__(self, fchk_path):\n",
    "        self.fchk_path = fchk_path\n",
    "        self.obj = fchk = FormchkInterface(fchk_path)\n",
    "        self.mol_weight = fchk.key_to_value(\"Real atomic weights\")\n",
    "        self.natm = natm = self.mol_weight.size\n",
    "        self.mol_coord = fchk.key_to_value(\"Current cartesian coordinates\").reshape((natm, 3))\n",
    "        self.mol_hess = fchk.hessian()\n",
    "        self.mol_hess = (self.mol_hess + self.mol_hess.T) / 2\n",
    "        self.mol_hess = self.mol_hess.reshape((natm, 3, natm, 3))\n",
    "        self._theta = NotImplemented     # Force constant tensor\n",
    "        self._proj_inv = NotImplemented  # Inverse space of translation and rotation of theta\n",
    "        self._freq = NotImplemented      # Frequency in cm-1 unit\n",
    "        self._q = NotImplemented         # Unnormalized normal coordinate\n",
    "        self._qnorm = NotImplemented     # Normalized normal coordinate\n",
    "        \n",
    "    @property\n",
    "    def theta(self):\n",
    "        if self._theta is NotImplemented:\n",
    "            natm, mol_hess, mol_weight = self.natm, self.mol_hess, self.mol_weight\n",
    "            self._theta = np.einsum(\"AtBs, A, B -> AtBs\", mol_hess, 1 / np.sqrt(mol_weight), 1 / np.sqrt(mol_weight)).reshape(3 * natm, 3 * natm)\n",
    "        return self._theta\n",
    "    \n",
    "    @property\n",
    "    def center_coord(self):\n",
    "        return (self.mol_coord * self.mol_weight[:, None]).sum(axis=0) / self.mol_weight.sum()\n",
    "    \n",
    "    @property\n",
    "    def centered_coord(self):\n",
    "        return self.mol_coord - self.center_coord\n",
    "    \n",
    "    @property\n",
    "    def rot_eig(self):\n",
    "        natm, centered_coord, mol_weight = self.natm, self.centered_coord, self.mol_weight\n",
    "        rot_tmp = np.zeros((natm, 3, 3))\n",
    "        rot_tmp[:, 0, 0] = centered_coord[:, 1]**2 + centered_coord[:, 2]**2\n",
    "        rot_tmp[:, 1, 1] = centered_coord[:, 2]**2 + centered_coord[:, 0]**2\n",
    "        rot_tmp[:, 2, 2] = centered_coord[:, 0]**2 + centered_coord[:, 1]**2\n",
    "        rot_tmp[:, 0, 1] = rot_tmp[:, 1, 0] = - centered_coord[:, 0] * centered_coord[:, 1]\n",
    "        rot_tmp[:, 1, 2] = rot_tmp[:, 2, 1] = - centered_coord[:, 1] * centered_coord[:, 2]\n",
    "        rot_tmp[:, 2, 0] = rot_tmp[:, 0, 2] = - centered_coord[:, 2] * centered_coord[:, 0]\n",
    "        rot_tmp = (rot_tmp * mol_weight[:, None, None]).sum(axis=0)\n",
    "        _, rot_eig = np.linalg.eigh(rot_tmp)\n",
    "        return rot_eig\n",
    "    \n",
    "    @property\n",
    "    def proj_scr(self):\n",
    "        natm, centered_coord, rot_eig, mol_weight = self.natm, self.centered_coord, self.rot_eig, self.mol_weight\n",
    "        rot_coord = np.einsum(\"At, ts, rw -> Asrw\", centered_coord, rot_eig, rot_eig)\n",
    "        proj_scr = np.zeros((natm, 3, 6))\n",
    "        proj_scr[:, (0, 1, 2), (0, 1, 2)] = 1\n",
    "        proj_scr[:, :, 3] = (rot_coord[:, 1, :, 2] - rot_coord[:, 2, :, 1])\n",
    "        proj_scr[:, :, 4] = (rot_coord[:, 2, :, 0] - rot_coord[:, 0, :, 2])\n",
    "        proj_scr[:, :, 5] = (rot_coord[:, 0, :, 1] - rot_coord[:, 1, :, 0])\n",
    "        proj_scr *= np.sqrt(mol_weight)[:, None, None]\n",
    "        proj_scr.shape = (-1, 6)\n",
    "        proj_scr /= np.linalg.norm(proj_scr, axis=0)\n",
    "        return proj_scr\n",
    "    \n",
    "    @property\n",
    "    def proj_inv(self):\n",
    "        if self._proj_inv is NotImplemented:\n",
    "            natm, proj_scr, theta = self.natm, self.proj_scr, self.theta\n",
    "            proj_inv = np.zeros((natm * 3, natm * 3))\n",
    "            proj_inv[:, :6] = proj_scr\n",
    "            x_col = 6 - 1\n",
    "            for A_col in range(6, natm * 3):\n",
    "                stat = True\n",
    "                while stat:  # first xcol - 6 values in vector t should be zero\n",
    "                    x_col += 1\n",
    "                    x0 = np.zeros((natm * 3, ))\n",
    "                    x0[x_col - 6] = 1\n",
    "                    t = x0 - proj_inv[:, :A_col].T @ x0 @ proj_inv[:, :A_col].T\n",
    "                    t /= np.linalg.norm(t)\n",
    "                    stat = (np.linalg.norm(t[:x_col - 6]) > 1e-7)\n",
    "                proj_inv[:, A_col] = t\n",
    "            proj_inv = proj_inv[:, 6:]\n",
    "            self._proj_inv = proj_inv\n",
    "        return self._proj_inv\n",
    "    \n",
    "    def _get_freq_qdiag(self):\n",
    "        natm, proj_inv, theta, mol_weight = self.natm, self.proj_inv, self.theta, self.mol_weight\n",
    "        e, q = np.linalg.eigh(proj_inv.T @ theta @ proj_inv)\n",
    "        freq = np.sqrt(np.abs(e * E_h * 1000 * N_A / a_0**2)) / (2 * np.pi * c_0 * 100) * ((e > 0) * 2 - 1)\n",
    "        self._freq = freq\n",
    "        q_unnormed = np.einsum(\"AtQ, A -> AtQ\", (proj_inv @ q).reshape(natm, 3, (proj_inv @ q).shape[-1]), 1 / np.sqrt(mol_weight))\n",
    "        q_unnormed = q_unnormed.reshape(-1, q_unnormed.shape[-1])\n",
    "        q_normed = q_unnormed / np.linalg.norm(q_unnormed, axis=0)\n",
    "        self._q = q_unnormed\n",
    "        self._qnorm = q_normed\n",
    "    \n",
    "    @property\n",
    "    def freq(self):\n",
    "        if self._freq is NotImplemented:\n",
    "            self._get_freq_qdiag()\n",
    "        return self._freq\n",
    "    \n",
    "    @property\n",
    "    def q(self):\n",
    "        if self._q is NotImplemented:\n",
    "            self._get_freq_qdiag()\n",
    "        return self._q\n",
    "    \n",
    "    @property\n",
    "    def qnorm(self):\n",
    "        if self._qnorm is NotImplemented:\n",
    "            self._get_freq_qdiag()\n",
    "        return self._qnorm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "fchk = FreqAnal(\"CH4.fchk\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "code_folding": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1675.76946, 1675.76946, 1675.76947, 1903.7231 , 1903.72311, 3525.8641 , 3786.57203, 3786.57206, 3786.57207])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fchk.freq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "31.22307017202692"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "FactIR = (e_c * N_A) * np.sqrt(1 / 10000000 * np.pi / 3)\n",
    "FactIR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.602176634e-19"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "e_c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 2.57742, -0.     , -0.     ],\n",
       "       [ 0.     ,  0.54676,  2.51876],\n",
       "       [ 0.     ,  2.51876, -0.54676],\n",
       "       [-0.     , -0.     ,  0.     ],\n",
       "       [-0.     ,  0.     , -0.     ],\n",
       "       [ 0.     ,  0.     ,  0.     ],\n",
       "       [ 0.     , -0.01668,  0.68369],\n",
       "       [ 0.68389, -0.     ,  0.     ],\n",
       "       [ 0.     , -0.68369, -0.01668]])"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dipderiv_mode = np.einsum(\"Ar, Aq -> qr\", fchk.obj.dipolederiv(), fchk.q) * FactIR\n",
    "dipderiv_mode"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "natm = 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[[-0.     , -0.     ,  0.     ],\n",
       "        [ 0.5    , -0.     ,  0.     ],\n",
       "        [-0.5    , -0.     ,  0.     ],\n",
       "        [-0.     ,  0.28868,  0.40825],\n",
       "        [ 0.     , -0.28868, -0.40825]],\n",
       "\n",
       "       [[ 0.     , -0.     ,  0.     ],\n",
       "        [-0.     , -0.5    ,  0.     ],\n",
       "        [ 0.     , -0.16667,  0.4714 ],\n",
       "        [-0.28868,  0.33333, -0.2357 ],\n",
       "        [ 0.28868,  0.33333, -0.2357 ]],\n",
       "\n",
       "       [[-0.     ,  0.     ,  0.     ],\n",
       "        [ 0.     , -0.     , -0.5    ],\n",
       "        [ 0.     ,  0.4714 ,  0.16667],\n",
       "        [-0.40825, -0.2357 ,  0.16667],\n",
       "        [ 0.40825, -0.2357 ,  0.16667]]])"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fchk.qnorm.reshape(natm, 3, 3 * natm - 6)[:, :, 3:6].transpose((2, 0, 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1675.76946, 1675.76946, 1675.76947, 1903.7231 , 1903.72311, 3525.8641 , 3786.57203, 3786.57206, 3786.57207])"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fchk.freq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "                     1                      2                      3\n",
    "                     A                      A                      A\n",
    " Frequencies --  1675.7695              1675.7695              1675.7695\n",
    " Red. masses --     1.1714                 1.1714                 1.1714\n",
    " Frc consts  --     1.9382                 1.9382                 1.9382\n",
    " IR Inten    --     6.6431                 6.6431                 6.6431\n",
    " Raman Activ --     3.8439                 3.8439                 3.8439\n",
    " Depolar (P) --     0.7500                 0.7500                 0.7500\n",
    " Depolar (U) --     0.8571                 0.8571                 0.8571\n",
    "  Atom  AN      X      Y      Z        X      Y      Z        X      Y      Z\n",
    "     1   6     0.12   0.00   0.00     0.00   0.02   0.12     0.00   0.12  -0.02\n",
    "     2   1    -0.60   0.00   0.00     0.00  -0.08   0.11     0.00  -0.60  -0.01\n",
    "     3   1    -0.60   0.00   0.00     0.00   0.23  -0.49     0.00   0.01   0.29\n",
    "     4   1    -0.12   0.28  -0.20    -0.16  -0.17  -0.53     0.30  -0.43  -0.05\n",
    "     5   1    -0.12  -0.28   0.20     0.16  -0.17  -0.53    -0.30  -0.43  -0.05"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([6.64311, 6.64311, 6.64311, 0.     , 0.     , 0.     , 0.46771, 0.46771, 0.46771])"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ir_intensities = (dipderiv_mode ** 2).sum(axis=1)\n",
    "ir_intensities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
